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ABSTRACT 

We present VULCAN/2D multi-group flux-limited-diffusion radiation hydrodynamics simulations of binary 
neutron star (BNS) mergers, using the Shen equation of state, covering >100ms, and starting from azimuthal- 
averaged 2D slices obtained from 3D SPH simulations of Rosswog & Price for 1 .4 M Q (baryonic) neutron stars 
with no initial spins, co-rotating spins, and counter-rotating spins. Snapshots are post-processed at 10 ms inter- 
vals with a multi-angle neutrino-transport solver. We find polar-enhanced neutrino luminosities, dominated by 
D e and neutrinos at peak, although v e emission may be stronger at late times. We obtain typical peak neu- 
trino energies for v e , v e , and "u^" of ~12, ~16, and ~22 MeV. The super-massive neutron star (SMNS) formed 
from the merger has a cooling timescale of <1 s. Charge-current neutrino reactions lead to the formation of a 
thermally-driven bipolar wind with <M>^10~ 3 M Q s" 1 , bary on-loading the polar regions, and preventing any 
production of a 7-ray burst prior to black-hole formation. The large budget of rotational free energy suggests 
magneto-rotational effects could produce a much greater polar mass loss. We estimate that < 10 -4 M Q of mate- 
rial with electron fraction in the range 0. 1-0.2 become unbound during this SMNS phase as a result of neutrino 
heating. We present a new formalism to compute the v,Vj annihilation rate based on moments of the neutrino 
specific intensity computed with our multi-angle solver. Cumulative annihilation rates, which decay as f~ l s , 
decrease over our 100 ms window from a few x 10 50 to ~10 49 erg s" 1 , equivalent to a few x 10 54 to ~10 53 e~e + 
pairs per second. 

Subject headings: stars: neutron - stars: supernovae: general - neutrinos - rotation - Gamma-ray: bursts - 
Hydrodynamics 



1. INTRODUCTION 

Coalescing neutron stars are one of the primary progenitor 
candidates for short-duration (i.e., < 2 s) 7-ray bursts (GRBs; 
Paczynski 1986; Goodman 1986; Eichleret al. 1989; Narayan 
et al. 1992; Mochkovitch et al. 1993). In-spiral of the two neu- 
tron star components occurs due to energy and angular mo- 
mentum loss through gravitational radiation (Taylor & Weis- 
berg 1982; Weisberg & Taylor 2005), which is emitted as the 
system evolves towards coalescence. Modern 7-ray and X- 
ray satellites have considerably improved our understanding 
of short-hard GRBs (Nakar 2007), and have more strongly as- 
sociated them with binary neutron star (BNS) merger events 
(Lee & Ramirez-Ruiz 2007). Gravitational wave detectors 
such as LIGO (Abramovici et al. 1992) will most likely im- 
prove our understanding of these events in the coming decade. 

The engine powering such short-hard GRBs and their as- 
sociated high-Lorentz-factor ejecta is thought to be linked to 
mass accretion from a torus onto a central black hole, formed 
when the super-massive neutron star (SMNS) born from the 
coalescence eventually collapses. In recent years, models 
of GRB production from a high-mass, magnetar-like, neu- 
tron star have also received some attention (Usov 1992, 1994; 
Kluzniak & Ruderman 1998; Rosswog et al. 2003; Dai et al. 
2006; Metzger et al. 2007; Dessart et al. 2007). We will ad- 
dress the likelyhood that the SMNS produces a GRB prior to 
black-hole formation in § 3.2. A fraction of the gravitational 
energy in the disk around this black hole is radiated as neutri- 
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nos which power a disk wind. A fraction of these neutrinos 
and antineutrinos annihilate into e~e + . Magnetic processes 
may, through their role in angular momentum transport, be 
decisive in setting the time scale between merger and collapse 
and, in addition, they may extract rotational energy from the 
central black hole (Blandford & Znajek 1977). The ejecta 
may be confined by the magnetic field morphology or by 
the intense neutrino-driven bary on-loaded wind thought to ac- 
company coalescence and black hole mass accretion (Levin- 
son & Eichler 2000; Rosswog & Ramirez-Ruiz 2003; Aloy 
et al. 2005). 

Numerical simulations of BNS mergers have a rich history, 
and have been considerably refined over the years, moving 
from 2D to 3D; from Newtonian to post-Newtonian, confor- 
mally flat, and finally to full General Relativity (GR); from 
simple polytropic to more sophisticated equations of state 
(EOS); from the neglect of neutrinos to approximate neutrino- 
trapping schemes. BNS merger studies yield predictions for 
1) their associated gravitational wave signals (Ruffert et al. 
1996; Calder & Wang 2002; Faber & Rasio 2002; Shibata 
& Uryu 2002; Shibata et al. 2003, 2005; Oechslin & Janka 
2007; Anderson et al. 2008b), 2) the production of short-hard 
GRBs (Ruffert et al. 1997; Ruffert & Janka 1999; Rosswog 
& Ramirez-Ruiz 2003; Rosswog et al. 2003; Rosswog 2005; 
Janka et al. 2006; Birkl et al. 2007), and 3) the pollution of 
the environment by r-process nuclei (Lattimer & Schramm 
1974, 1976; Symbalisty & Schramm 1982; Eichler etal. 1989; 
Davies et al. 1994; Freiburghaus et al. 1999; Rosswog et al. 
1999; Surman etal. 2008). 

Different modeling ingredients and approaches have been 
employed. 3D Newtonian simulations without neutrino trans- 
port were performed by Oohara & Nakamura (1989, 1990), 
Nakamura & Oohara (1989, 1991), Zhuge et al. (1994, 1996), 
and Rosswog et al. (2000). Ruffert et al. (1996, 1997) 
included neutrinos using a grid-based code, and this was 
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followed by Rosswog & Davies (2002) and Rosswog & 
Liebendorfer (2003) using a 3D SPH code. Recently, Se- 
tiawan et al. (2004, 2006) applied such a neutrino-leakage 
scheme in their study of torus-disk mass accretion around a 
black hole formed in a BNS merger, and addressed neutrino 
emission and annihilation. Magnetic fields in BNS merger 
evolution were introduced in 3D Newtonian simulations by 
Price & Rosswog (2006) and Rosswog & Price (2007). 

In parallel, there have been improvements to the New- 
tonian approach to incorporate the effects of GR (although 
largely neglecting microphysics), which become important 
as the two neutron stars come closer and eventually merge. 
Post-Newtonian simulations were performed by Oohara & 
Nakamura (1992), Faber & Rasio (2000, 2002), Faber et al. 

(2001) , Ayal et al. (2001), and Calder & Wang (2002). 
The conformally-flat approximation was introduced in Wil- 
son et al. (1996) and followed by Oechslin et al. (2002), Faber 
et al. (2004), Oechslin & Janka (2006, 2007), and Oechslin 
et al. (2007). Full 2D GR simulations, sometimes including 
magnetic fields, have been performed for super-massive neu- 
tron stars (SMNS), possibly resulting from BNS mergers, by 
Baumgarte et al. (2000), Morrison et al. (2004), Duez et al. 
(2004, 2006), and Shibata et al. (2006), while full 3D GR 
simulations of the merger were carried out by Shibata & Uryu 

(2002) ; Shibata et al. (2003, 2005); Shibata & Uryu (2007), 
Anderson et al. (2008a), and Baiotti et al. (2008). 

All these simulations have been conducted at various levels 
of sophistication for the thermodynamic properties of matter, 
ranging from simplistic and not so physically-consistent poly- 
tropic EOSs, to those employing a detailed microphysical rep- 
resentation of nuclear matter at an arbitrary temperature (Lat- 
timer&Swesty 1991; Shen et al. 1998a,b). Investigations per- 
formed with such detailed EOSs have been conducted by, e.g., 
Ruffert et al. (1996, 1997), Rosswog et al. (1999), and Ross- 
wog & Davies (2002), and the dependency of BNS merger 
properties on the adopted EOS has been discussed by Oech- 
slin & Janka (2007) and Oechslin et al. (2007). Variations 
of a few times 10% in the maximum allowed mass are seen, 
depending on the compressiblity of nuclear matter, but re- 
cent observations may suggest that neutron stars with a grav- 
itational mass in excess of 2M Q do exist (Nice et al. 2005; 
Freire et al. 2007), supporting a stiff EOS for nuclear matter, 
such as the Shen EOS we employ here. In Fig. 1, we show, 
as a function of the central density, the baryonic and gravita- 
tional neutron star masses that obtain for our implementation 
of the Shen EOS. The maximum occurs at^l.3xl0 15 g cm" 3 , 
corresponding to a gravitational (baryonic) mass of 2.32 M Q 
(2.77 M Q ). 5 Moreover, in the context of differentially rotat- 
ing (and possibly magnetized) SMNSs, the maximum mass 
that can be supported by a given EOS may be increased by 
up to ^50% compared to the equivalent non-rotating object 
(Baumgarte et al. 2000; Morrison et al. 2004; Duez et al. 2004, 
2006; Shibata et al. 2006; Sumiyoshi et al. 2007; Kiuchi & 
Kotake 2008). Ultimately, understanding the mechanisms and 
timescales for angular momentum to be redistributed to lead 
to solid-body rotation is, therefore, important. One possible 
agency of redistribution is the magneto-rotational instability 
(Balbus & Hawley 1991a; Akiyama et al. 2003; Pessah et al. 
2006). The implications are non-trivial because the efficiency 
of angular-momentum transport determines in part whether a 
black hole forms promptly, or after a short or a long delay. It 
also determines how much mass can be accreted, i.e., whether 

5 Note that the values given in Dessart et al. 2008 are not exact. 
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FIG. 1. — Gravitational (solid line) and baryonic (dotted line) masses, 
as well as the corresponding neutron star radius (dashed line), versus cen- 
tral density derived for our Shen EOS and the Oppenheimer-Volkov equa- 
tion (one form of the general relativistic equation of hydrostatic equilibrium), 
and assuming a non-rotating neutron star characterized by a uniform elec- 
tron fraction of 0. 1 and a uniform temperature of 1 MeV. Here, the transi- 
tion to a black hole will occur when the central density in the neutron star 
reaches ~1.3 X 10 15 g cm -3 , corresponding to a gravitational (baryonic) mass 
of 2.32 Mq (2.77 Mq ) and a neutron star radius of ~ 1 3 km. 

there is ^0.01 or ^0.1 Mq available in the torus surrounding 
the black hole after it has formed, and what the timescale is 
over which such accretion can take place to power relativistic 
ejecta. Although indirect, the relevance to short-hard GRBs 
and their properties is obvious, and perhaps central. In fact, 
this issue is also germane to the production of collapsars and 
long-duration GRBs (Dessart et al. 2008). 

In this work, we perform 2D multi-group, flux-limited- 
diffusion (MGFLD), radiation hydrodynamics of merged 
BNSs using the Shen EOS. The main goal of this work is to 
document in detail the neutrino signatures from such merger 
events. Our approach is to solve the neutrino transport prob- 
lem self-consistently (although with a flux-limiter and assum- 
ing diffusion), in combination with the dynamics of the sys- 
tem (but assuming axisymmetry), over >100ms. At selected 
times, we post-process such models with the multi-angle, S n , 
neutrino transport algorithm of Livne et al. (2004) and de- 
scribed recently in Ott et al. (2008). In particular, our in- 
vestigation allows for the first time the computation of the 
neutrino-antineutrino annihilation rate from a solution of the 
transport equation, rather than using an approximate leakage 
scheme. Our investigation applies strictly to the neutron-star 
phase of such BNS mergers. 

This paper is organized as follows. In § 2, we present the 
initial models we employ in our work, based on 3D SPH sim- 
ulations of BNS mergers with different spin configurations 
(§ 2.1). In § 2.2, we present the characteristics and proce- 
dures we employ to evolve such initial conditions with the 
MGFLD radiation hydrodynamics code VULCAN/2D. In § 3, 
we present our results for the neutrino signatures (§ 3.1) and 
the dynamics of BNS mergers (§ 3.2). In § 4, we address 
the neutrino-antineutrino annihilation rate in our three BNS 
merger models. We first present results employing the ap- 
proach used so far (Ruffert et al. 1996, 1997; Rosswog & 
Liebendorfer 2003; Setiawan et al. 2004, 2006; Birkl et al. 
2007), based on a leakage scheme for the neutrino-flavor de- 
pendent opacities and emissivities and a summation of all 
paired grid cells contributing at any given location. In § 4.2, 
we then present a new formalism which uses moments of the 
neutrino specific intensity computed with a multi-angle, S n , 
scheme. Finally, we present our conclusions in § 5 and dis- 
cuss the most striking implications for the powering of short- 
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duration GRBs. For completeness, in the appendix, we apply 
our formalism to the computation of the annihilation rate in 
the context of slow- and fast-rotating single protoneutron stars 
(PNSs). 

2. MODEL DESCRIPTION 

The main goal of this work is to understand the long- 
term (over many rotation periods), post-coalescence, evolu- 
tion of BNS mergers, with particular attention to neutrino 
signatures (flavor dependence, angular distribution, and ra- 
dial distribution). Starting from azimuthal-averaged slices 
constructed from 3D Smooth-Particle-Hydrodynamics (SPH) 
simulations with the MAGMA code (Rosswog & Price 
2007) and taken a few milliseconds after coalescence, we 
perform two-dimensional (axisymmetric) multi-group flux- 
limited-diffusion (MGFLD) radiation hydrodynamics simula- 
tions using the code VULCAN/2D (Livne et al. 2004; Dessart 
et al. 2006a; Burrows et al. 2007b; Ott et al. 2008). Below, we 
present the properties of the three BNS merger configurations 
from which we start (§ 2.1), and then describe our approach 
with VULCAN/2D in more detail (§ 2.2). 

2. 1 . Initial conditions 

The MAGMA code (Rosswog & Price 2007) is a state-of- 
the-art SPH code that contains physics modules that are rel- 
evant to compact binary mergers and on top implements a 
slew of numerical improvements over most "standard" SPH 
schemes. A detailed code-description can be found in Ross- 
wog & Price (2007), while results are presented in Price & 
Rosswog (2006) and Rosswog (2007). 

We briefly summarize the most important physics modules. 
For the thermodynamic properties of neutron star matter we 
use a temperature-dependent relativistic mean-field equation 
of state (Shen et al. 1998a,b). It can handle temperatures 
from to 100 MeV, electron fractions from T e = (pure neu- 
tron matter) up to 0.56 and densities from about 10 to more 
than 10 15 gem" 3 . No attempt is made to include matter con- 
stituents that are more exotic than neutrons and protons at 
high densities. For more details we refer the reader to Ross- 
wog & Davies (2002). 

The MAGMA code contains a detailed multi-flavor neu- 
trino leakage scheme. An additional mesh is used to calculate 
the neutrino opacities that are needed for the neutrino emis- 
sion rates at each particle position. The neutrino emission 
rates are used to account for the local cooling and the com- 
positional changes due to weak interactions such as electron 
captures. A detailed description of the neutrino treatment can 
be found in Rosswog & Liebendorfer (2003). 

The self-gravity of the fluid is treated in a Newtonian fash- 
ion. Both the gravitational forces and the search for the par- 
ticle neighbors are performed with a binary tree that is based 
on the one described in Benz et al. (1990). These tasks are the 
computationally most expensive part of the simulations and 
in practice they completely dominate the CPU-time usage. 
Forces emerging from the emission of gravitational waves are 
treated in a simple approximation. For more details, we re- 
fer to the literature (Rosswog et al. 2000; Rosswog & Davies 
2002) 

In terms of numerical improvements over "standard" SPH 
techniques, the code contains the following: 

• To restrict shocks to a numerically resolvable width ar- 
tificial viscosity is used. The form of the artificial vis- 
cosity tensor is oriented at Riemann solvers (Monaghan 



1997). The resulting equations are similar to those con- 
structed for Riemann solutions of compressible fluid 
dynamics. In order to apply the artificial viscosity terms 
only where they are really needed, i.e., close to a shock, 
the numerical parameter that controls the strength of the 
dissipative terms is made time dependent, as suggested 
by Morris & Monaghan (1997). An extra equation is 
solved for this parameter which contains a source term 
that triggers on the shock and a term causing an expo- 
nential decay of the parameter in the absence of shocks. 
Tests can be found in Rosswog et al. (2000) and an il- 
lustration of the time-dependent viscosity parameter in 
the context of Sod's shock tube problem can be found 
in Rosswog et al. (2008, their Fig. 1). 

• A self-consistent treatment of extra terms in the hydro- 
dynamics equations that arise from varying smoothing 
lengths (so-called "grad-h"-terms; Springel & Hern- 
quist 2002; Monaghan 2002; Price 2004; Rosswog & 
Price 2007). 

• A consistent implementation of adaptive gravitational 
softening lengths as described in Price & Monaghan 
(2007). 

• The option to evolve magnetic fields via so-called Eu- 
ler potentials (Stern 1970) that guarantee that the con- 
straint V • B = is fulfilled. The details can be found in 
Rosswog & Price (2007) and Rosswog & Price (2008). 

Our work starts from 3D SPH simulations of BNS mergers 
in an isolated binary system, made up of two individual neu- 
tron stars, each with a baryonic mass of 1.4M Q (equivalent 
to a gravitational mass of ~1.3 M Q ), and separated by a dis- 
tance of 48 km. Their initial hydrostatic structure is computed 
by solving the Newtonian equations. The properties of nu- 
clear matter (pressure, energy, entropy etc.) are obtained by 
interpolation in the Shen EOS (Shen et al. 1998a,b), assuming 
zero temperature and (3 equilibrium. A large number of SPH 
particles (N f=a 10 6 ) are then mapped onto the resulting den- 
sity profiles These particles are further relaxed to find their 
true equilibrium state (see, e.g., Rosswog & Price 2007). 

Because the tidal synchronization time exceeds the grav- 
itational decay time, BNSs cannot be tidally locked during 
their inspiral phase, and are thus expected to be close to irro- 
tational (see, e.g. Bildsten & Cutler 1992 and Kochanek 1992. 
However, past investigations of the coalescence of BNSs sug- 
gest that the intrinsic spin of each neutron star is an impor- 
tant parameter, as it determines, for example, the width of the 
baryon-poor funnel above the forming central object that is 
thought to play a decisive role in the launching of a GRB (see, 
e.g. Zhugeetal. 1996;Faberetal. 2001; Oechslin et al. 2007). 
Therefore, to encompass the range of possible outcomes, we 
simulate three initial configurations: 1) no spins (no intrin- 
sic neutron star spin), 2) co-rotating spins (each neutron star 
has a period equal to that of the orbit, and rotates in the same 
direction), and 3) counter-rotating spins (same spin period as 
for co-rotating spins, but with rotation in a reverse direction 
to that of the orbit; the spins are counter-aligned). In the rest 
of this paper, we broadly refer to these models as the no-spin, 
the co-rotating spin, and the counter-rotating spin BNS mod- 
els. The co-rotating and counter-rotating cases are included 
here to bracket the range of possibilities, but they represent 
unlikely extremes. We leave to future work the considera- 
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FIG. 2. — Colormap of the temperature (MeV; log scale) at (top row; inner 150x 150km 2 region) and 100 ms (bottom row, inner 250x250 km 2 region) after 
the start of the simulation. In the bottom panels, velocity vectors are added, with a saturation length at 15% the width of the display and corresponding to a 
velocity of 30,000 km s . From left to right columns, we show the no-spin, co-rotating, and counter-rotating models. We also overplot white density contours 
for every decade, starting at a maximum of 10 14 gcirT 3 . Note that the minimum of the colorbar changes between the top and the bottom rows. 

MGFLD radiation-hydrodynamics code VULCAN/2D (Livne 
et al. 2004; Dessart et al. 2006a; Burrows et al. 2007b). The 
gravitational potential is assumed Newtonian and computed 
using a multipole solver (Dessart et al. 2006a). Due to the 
difficulty of handling the superluminal Alfven speeds prevail- 
ing in the low-density material surrounding the BNS merger 
for even modest field strengths, we do not investigate the de- 
pendence on magnetic field strength and morphology. How- 
ever, our study presents a self-consistent treatment of neu- 
trino emission, scattering, and absorption in a multi-species 
MGFLD context. Associated neutrino-matter coupling source 
terms are included in the momentum and energy conserva- 
tion equations and their relevance to the merger evolution is, 
thus, computed explicitly, for a typical duration of > 100 ms. 
For these dynamical calculations, we include all neutrino pro- 
cesses described in Burrows et al. (2006), but neglect the sec- 
ondary processes of neutrino-electron scattering. The trans- 
port solution in the MGFLD scheme we employ is solved for 
the three neutrino species v e , v e , and (which groups to- 
gether the fj, and r neutrinos) and at eight neutrino energies: 
2.50, 6.87, 12.02, 21.01, 36.74, 64.25, 112.36, 196.48 MeV. 
The energy spacing is constant in the log. 

The choice of spatial grid is determined primarily by the 
very aspherical density distribution of the merger and the very 
steep density drop-off at the neutron star surface in the po- 
lar direction. We present in Figs. 2-3 (top row) the initial 
temperature (logarithmic scale and in MeV), electron frac- 
tion, and the density distribution (white line contours over- 
plotted for every decade between 10 7 and 10 14 gcm" 3 ) for 
each BNS merger configuration that we map onto the VUL- 
CAN/2D grid. Note that fast rotation in the inner region, 



tion of higher mass BNS mergers (Rosswog & Ramirez-Ruiz 
2003) or asymmetric systems (Rosswog et al. 2000). 

In practice, we post-process these 3D SPH simulations 
by constructing azimuthal-averages, selecting a time after 
coalescence (when the two neutron stars come into con- 
tact) of 10.0 (no-spins), 12.4 (co-rotating spins), and 20.4 ms 
(counter-rotating spins), corresponding in each VULCAN/2D 
simulations to the time origin. At these times, these systems 
are evolving towards, but have not yet reached, complete ax- 
isymmetric configurations. This is a compromise made to 
model with VULCAN/2D the epoch of peak neutrino bright- 
ness. At these initial times, the amount of mass with a density 
greater than 10 13 g cm" 3 is 2.5 M in the no-spin BNS model, 
2.31 M Q in the co-rotating spin model, and 2.43 M Q in the 
counter-rotating spin model. In the same model order, and 
at these times, the mass contained inside the ~13km radius 
of the highest mass neutron star allowed by the Shen EOS 
is only 1.11, 1.26, and 1.03 M Q . The SPH simulations em- 
ployed to generate these models are Newtonian, and, thus, 
one would expect more mass at small radii and systematically 
higher densities if allowance for GR effects were made, but 
the use of the highly-incompressible Shen EOS together with 
the significant amount of mass in quasi-Keplerian motion and 
on wide orbits suggests that the systems should survive at least 
for some time before the general-relativistic gravitational in- 
stability occurs. 

2.2. VULCAN/2D 

Starting from fluid variables constructed with azimuthal av- 
erages of the 3D SPH simulations described above, we follow 
the evolution of three BNS merger configurations with the 
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FIG. 3. — Same as for Fig. 2, but now for the electron fraction F c . Note that the extent of the displayed region and the extrema of the colorbar change between 
the top and the bottom rows. In the counter-rotating spin model, an axis problem develops early on in the VULCAN/2D simulation (visible as a spurious low Y e 
"needle" along the rotation axis), despite the high spatial resolution employed. 



although sub-Keplerian, displaces the density maximum by 
8 km from the center and along the equator in the no-spin BNS 
model. These density peaks also correspond to extrema in the 
electron fraction of ~0.1 at this time. To setup the hybrid 
VULCAN/2D grid, we position the transition radius between 
the inner Cartesian and the outer spherical polar regions at 
12 km. The density at this location is at all times in excess of 
10 14 g cm" 3 , and thus - due to the stiffness of the Shen EOS in 
that regime - away from the regions of large density gradient. 
The resolution in this inner region is typically 200 x200m 2 , 
corresponding to 12/0.2 = 60 zones in each (cylindrical) di- 
rection r and z from the center to the transition radius. Be- 
yond the transition radius, we use a logarithmically increas- 
ing radial grid spacing with Ar/r = 1.85%, using 301 zones 
to 3000 km. The grid covers one hemisphere, from the rota- 
tion axis to the equator (both treated as reflecting boundaries), 
and the computation assumes axisymmetry, i.e., the azimuthal 
gradients are zero. We fill the grid outside the merger with 
material having a low density of 5 x 10 3 g cm" 3 and a low tem- 
perature of 4x 10 8 K. Below Y e =0.05, we compute thermody- 
namic variables, as well as opacities/emissivities, by adopting 
an electron fraction of 0.05. Given the smooth variation for 
the corresponding quantitites (pressure, entropy, mean-free 
path, etc.) at this level of neutron richness, this approxima- 
tion is expected to be quite good. 

After the MGFLD radiation-hydrodynamics simulations 
are completed, we post-process individual snapshots keep- 
ing the hydrodynamics frozen, and relaxing only the radia- 
tion quantities. This operation is performed for each model 
at 10 ms intervals over the whole evolution using the multi- 
angle, S„, radiation-transport module discussed in Livne et al. 
(2004) and described more recently in Ott et al. (2008). We 



refer the reader to those papers for details on the method and 
the nomenclature. This S„ variant yields the explicit angular 
dependence of the neutrino specific intensity, and, thus, rep- 
resents a more accurate modeling of the strongly anisotropic 
neutrino luminosity in such highly aspherical systems (Ott 
et al. 2008). The S„ calculation is done with the same num- 
ber of energy groups (i.e., eight) and on the same spatial grid. 
In the S„ algorithm, the transition at depth to MGFLD de- 
scribed in Ott et al. (2008) is natural here, since it occurs 
at 12 km, and, thus, in regions where the densities are nu- 
clear. However, the spatial resolution is not fully satisfac- 
tory at the neutron star surface and along the polar direc- 
tion, a region where the density decreases by few orders of 
magnitude in just a few kilometers. This, thus, represents 
a challenge for radiation transport on our Eulerian grid. In 
the first instance, we relax the MGFLD results with S„ adopt- 
ing 8 ^-angles. In the S„ approach, n such ^-angles translate 
into «(« + 2)/2 directions mapping the unit sphere uniformly. 
Once converged, we perform an accuracy check by remap- 
ping the angle-dependent neutrino radiation field from 8 to 16 
^-angles, and then relaxing this higher angular resolution sim- 
ulation. Due to the increased computational costs, we use the 
high (16 ^-angles) angular-resolution configuration only for 
snapshots at 10, 60, and 100 ms, but for all three BNS merger 
configurations studied here. Importantly, the explicit angle- 
dependence of the neutrino-radiation field allows us to com- 
pute its various moments, which we employ in a novel formu- 
lation of the neutrino-antineutrino annihilation rate (§ 4.2). 
We find that these results compare favorably with the simpler 
approach of Ruffert et al. (1996, 1997), although the multi- 
angle simulations offer a number of interesting new insights, 
which we explore in § 4. 
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FIG. 4. — Time evolution of the density distribution along the equatorial (top row) and the polar (bottom row) directions for the BNS merger models with 
initially no spins (left), co-rotating spins middle), and counter-rotating spins (right). 



3. BNS MERGER EVOLUTION AND NEUTRINO SIGNATURES 

In each model, the restart with VULCAN/2D is followed by 
a transient phase that lasts a few milliseconds and that is char- 
acterized by high-frequency oscillations of the highest density 
material, located within 10-20 km of the center. Given the 
good match between the Shen EOS used by Rosswog and in 
this work (pressures differ by at most a few percent), we spec- 
ulate that these oscillations are caused in part by the glitch 
introduced through the azimuthal averaging of the 3D SPH 
simulation snaphsot. However, they may also reflect a funda- 
mental dynamical property of this early phase. Similar oscil- 
lations are seen in 3D simulations performed with MAGMA, 
and have also been reported in the literature, e.g. by Baiotti 
et al. (2008). The associated shocks generate thermal energy 
(note that the thermal part of the pressure is sub-dominant 
at nuclear densities), but this thermal component is negligible 
compared, for example, with what is needed to power the neu- 
trino luminosities we see. Any initial mismatch is thus merely 
a small transient which has a negligible impact on the long- 
term evolution of the BNS mergers we study in this work. 

The evolutions of the BNS merger models with initially no 
spins, co-rotating spins, and counter-rotating spins, are quali- 
tatively similar. As discussed above, at the start of the VUL- 
CAN simulations, the matter distribution is far from station- 
ary, with a significant amount of material at sub-nuclear den- 
sities and at large distances from the center. Throughout the 
~ 100 ms of evolution we follow, matter in the inner few hun- 
dred kilometers settles in and comes to rest at the neutron star 
surface, while material beyond a few hundred kilometers and 
located at low, near equatorial latitudes, migrates outward due 
to the large thermal pressure gradient and the strong centrifu- 
gal effects. This expansion also takes place in the vertical di- 
rection, leading to the formation of a low-density cocoon sur- 
rounding an inner and denser disk. In the intermediate region, 
i.e., at radii of 50-100 km and along the equatorial direction, 
material evolves towards quasi-Keplerian motion, with little 
inward or outward migration. In Figs. 2-3, we show the T- 
Y e distributions at the start of the VULCAN/2D simulations 
(top row), and 100 ms later (bottom row), together with iso- 
density contours (shown in white) for every decade starting at 
a maximum value of 10 14 g cm" 3 . 

At 100 ms after the start of the VULCAN/2D simulations, 



each merger has reached a quasi-steady equilibrium configu- 
ration, characterized by a large equator to pole radius ratio of 
5/1, visible from the latitudinal variation of the extent of the 
10 10 gcm" 3 density contour (bottom row panels in Fig. 2-3) 
or from radial slices of the density in the polar and equatorial 
directions (Fig. 4). This is a general result for fast-rotating 
neutron stars (or other degenerate objects like white dwarfs) 
which has been documented in various contexts in the past 
(Ostriker & Mark 1 968 ; Hachisu 1 986; Liu & Lindblom 200 1 ; 
Yoon & Langer 2005; Rosswog & Davies 2002; Kiuchi & Ko- 
take 2008). At the start of the VULCAN simulations, high- 
density material was located in a thin (<20km) structure ex- 
tending < 300 km in the equatorial direction. After 100 ms, 
90% of the total mass is contained within 30 km of the cen- 
ter (corresponding to ^2.5 M Q ) and at densities greater than 
10 13 gcm" 3 . At the same time, the remaining ~10% (corre- 
sponding to ~0.2M Q ) is located in a quasi-Keplerian disk 
with an outer edge at <100km (Fig. 4-5; see also Fig. 12). 
This is also illustrated in Fig. 6 for the no-spin model and at 
60 ms after the start of the VULCAN/2D simulation, a figure 
in which we show the interior baryonic mass as a function of 
spherical radius, and how it compares with various multiples 
of the Schwarschild radius. 

3.1. Neutrino signatures 

Neutrino processes of emission and absorption/scattering 
in BNS merger simulations have been treated using leak- 
age schemes by Ruffert et al. (1996, 1997), Rosswog & 
Liebendorfer (2003), Setiawan et al. (2004, 2006) and Birkl 
et al. (2007). In this approach, the difficult solution of 
the multi-angle multi-group Bolzmann transport equation is 
avoided by introducing a neutrino loss timescale, for both en- 
ergy and number, and applied to optically-thick and optically- 
thin conditions. In addition, an estimated rate of electron-type 
neutrino loss allows the electron fraction to be updated. Over- 
all, benchmarking of these leakage schemes using more so- 
phisticated ID core-collapse simulations that solve the Boltz- 
mann equation, suggests that the resulting neutrino luminosi- 
ties are within a factor of a few at most of what would be 
obtained using a more accurate treatment. 

The work presented here offers a direct test of this. More- 
over, because the neutrino source terms are included in the 
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FIG. 5. — Time evolution of the angular velocity distribution along the equatorial direction for the BNS merger models with initially no spins (left), co-rotating 
spins middle), and counter-rotating spins (right). The dashed line gives the local Keplerian angular velocity (adopting the mass of the spherical volume interior to 
a given radius), suggesting that the material at low latitudes and located between ~30 and ~ 100 km forms a quasi-Keplerian disk, with densities on the order of 
10 12 gcm~ 3 . Note also the strong degree of differential rotation in the no-spin and counter-rotating spin cases, which is preserved throughout the VULCAN/2D 
simulation. 
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FIG. 6. — Interior baryonic mass versus spherical radius for the BNS 
merger model with no initial spins. The time corresponding to these data 
is 60 ms after the start of the VULCAN/2D simulation. Following Ruffert 
et al. (1996), we overplot the Schwarzschild radius R s (m(r)), scaled by 0.5, 
2.5, and 3, corresponding to the range of masses plotted on the vertical axis. 
Our results differ from those of Ruffert et al. (1996) in that only a small mass 
range falls within the radius 3xS s (m(r)) of the last stable circular orbit for 
an equal mass binary in harmonic coordinates (Kidder et al. 1992; Wex & 
Schaefer 1993). This difference likely results from the smaller mass binary 
system (2x 1.4Mq here compared to 2x 1.65 Mq in their work; the mass is 
the baryonic mass in both cases here) and the stiffer EOS we employ (Shen 
EOS versus Lattimer & Swesty EOS), with an incompressiblity of 180MeV. 

momentum and energy equations solved by VULCAN/2D, 
we can model the birth and subsequent evolution of the re- 
sulting neutrino-driven wind. Our simulations being 2D and 
less CPU intensive than 3D SPH simulations, can also be ex- 
tended to > 100 ms, and, thus, we can investigate the long-term 
evolution of the merger. In particular, we directly address the 
evolution of the neutrino luminosity and confront this with the 
results using the expedient of a steady-state often employed in 
work that is also limited to a short time span of 10-20ms after 
the merger. 

In the leakage schemes of Ruffert et al. (1996, 1997) and 
Rosswog & Liebendtirfer (2003), the neutrino emission pro- 
cesses treated are electron and positron capture on protons 
(yielding v e ) and neutrons (yielding v e ), respectively, and 
electron-positron pair annihilation and plasmon decay (each 
yielding electron-, /i-, and r-type neutrinos). For electron- 
type neutrinos, the dominant emission processes are the 
charged-current /^-processes. For neutrino opacity, Ruffert 
et al. include neutrino scattering on nucleons. Rosswog & 
Liebendorfer (2003) also treat electron-type neutrino absorp- 
tion on nucleons. 

The neutrino emission, absorption, and scattering pro- 
cesses used in VULCAN/2D simulations are those summa- 
rized in Burrows et al. (2006), and include all the above, plus 
electron neutrino absorption on nuclei and nucleon-nucleon 
bremsstrahlung. The latter is the dominant emission process 



of and v T neutrinos in PNSs (Thompson et al. 2000). In 
VULCAN/2D, the and v T neutrinos are grouped together 
and referred to as "i/^" neutrinos. We present in Fig. 7 the 
evolution of the neutrino luminosity for the BNS merger mod- 
els with initially no spins (top), co-rotating spins (center), 
and counter-rotating spins (bottom) over a typical timespan of 
> 1 00 ms, and for an adopted radius of 200 km along the equa- 
torial (polar) directions in black (red). In each case, we plot 
the v e (solid line), the D e (dashed line), and the "u^" (dash- 
dotted line) neutrinos. For comparison with the predictions 
of the leakage scheme of Ruffert et al. (1996), we have im- 
plemented their formalism, following exactly the description 
given in their paper in their Appendix A & B. We plot the cor- 
responding results at 10 ms intervals for each case (diamonds: 
v e ; triangles: v e ; squares: "v^"; see also § 4.1). 

For the BNS merger models with initially no spins, co- 
rotating spins, and counter-rotating spins, the neutrino sig- 
nal predicted by VULCAN/2D follows a qualitatively simi- 
lar evolution. The initial fast rise represents the time it takes 
the neutrino emission processes to ramp up to their equi- 
librium rates, modulated by the diffusion time out of the 
opaque core (and out of the not-so-opaque surface layers) 
and the light travel time of ~lms to the radius of 200 km 
where the luminosity is recorded. The peak luminosity oc- 
curs at about 5 ms after the start of the VULCAN/2D simu- 
lations, and is dominated by v e and "v^" neutrinos. The v e 
neutrino luminosity is systematically sub-dominant initially 
compared with the other two neutrino species. This is pri- 
marily a result of the high neutron richness of the merged 
object. Among the models, the stronger neutrino signal, in 
particular for the "v^' neutrinos, is produced in the counter- 
rotating model, which achieves the largest central tempera- 
ture (~ 17 MeV), followed by the no-spin model (~15 MeV). 
The co-rotating spin model has the weakest neutrino sig- 
nal of all three, mainly because the tidal locking leads to 
a very smooth merger and correspondingly low tempera- 
tures (the central temperature is about 10 MeV). Specifically, 
the angle-integrated, species-integrated, peak neutrino lumi- 
nosity for the BNS model with no initial spin is 1.5 xlO 53 
(2.2x 10 52 , 7.0x 10 52 , and 6.7x 10 52 ), for the co-rotating spin 
model 1.2x 10 53 (2.2 x 10 52 , 5.7 x 10 52 , and 3.9 x 10 52 ), and for 
the counter-rotating spin model 2.2 x 10 53 erg s" 1 (2.9xl0 52 , 
8.2xl0 52 , and 1.2xl0 53 ). The values given in parentheses 
correspond in each case to the angle-integrated peak lumi- 
nosity for the v e , the v e , and the "ly," neutrinos, respec- 
tively. Apart from a significant discrepancy with the "z^"- 
neutrino luminosities, the peak neutrino luminosities using 
our MGFLD approach compare well with those published in 
the literature based on a leakage scheme (Ruffert et al. 1997; 
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FIG. 7. — Time evolution of the MGFLD energy-integrated fluxes at 
R =200 km, scaled by a factor 4nR~ (thus equivalent to a luminosity), for 
the u e (solid), the v e (dashed), and the u Va" (dash-dotted) neutrinos, plot- 
ted along the equatorial (black) and polar (red) directions, and for the BNS 
merger models initially with no spins (top), co-rotating spins (middle), and 
counter-rotating spins (bottom). We also overplot the neutrino luminosities 
(diamonds: v e ; triangles: v e ; squares: 'V M ") computed in § 4.1 using the 
leakage scheme of Ruffert et al. (1996), but only for snapshots at 10ms in- 
tervals. Nucleon-nucleon bremsstrahlung processes are not treated in their 
scheme and this is the likely cause of the very low 'V„" neutrino luminos- 
ity compared to the VULCAN/2D prediction. Note that the fast rise of the 
VULCAN/2D neutrino luminosities reflects in part the light-travel time of 
~ 1 ms to the 200 km radius where the luminosity is recorded. Note that the 
luminosity unit used is the Bethe, i.e. 10 51 erg = 1 Bethe [1 B]. 



FIG. 8.— MGFLD flux spectra at R =200 km, scaled by a factor 4nR 2 
(thus equivalent to a luminosity spectrum), for v e (solid), v e (dashed), and 
"Uft" (dash-dotted) neutrinos, for the BNS merger models initially with no 
spins (left), co-rotating spins (middle), and counter-rotating spins (right), and 
plotted along the equatorial (black) and polar (red) directions. For all models, 
the time corresponds to 60 ms after the start of the simulation. 

Rosswog & Liebendorfer 2003). However, in all cases, cool- 
ing of the BNS merger causes a decrease of all neutrino lumi- 
nosities after the peak (i.e., rather than a long-term plateau), 
with a faster decline rate in the no-spin and counter-rotating 
spin cases. The luminosity decay rate is faster for the "i^" 
neutrinos than for the electron-type neutrinos. This translates 
into a strong decrease of the neutrino-antineutrino annihila- 
tion rate due to its scaling with the square of the neutrino 
luminosity (see § 4). In agreement with past work that fo- 
cused on the initial 10-20ms after the onset of coalescence, 
we find that the electron antineutrino luminosity dominates 
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FIG. 9. — Angular variation from pole to equator of the energy-integrated 
neutrino fluxes at 100 km (this radius is chosen smaller than for Fig. 7 to 
better reveal the flux anisotropy), scaled by a factor 4nR- (thus equivalent 
to a luminosity), for the v e (solid), the v e (dashed), and the 'V M " (dash- 
dotted) neutrinos, for the BNS merger models initially with no spins (top), co- 
rotating spins (middle), and counter-rotating spins (bottom). For each model, 
we show the MGFLD (black) and S n (red; using the 16 i?-angle calculations) 
predictions at 60 ms after the start of each simulation. Note the S„ oscillatory 
features in the neutrino flux, more prominant along the polar direction, which 
is an artifact of the S„ scheme. The range of the ordinate axis varies between 
frames. 



over the electron neutrino luminosity, but in contrast with the 
work of Ruffert et al. (1996) and Rosswog & Liebendorfer 
(2003), our "Vft" neutrino luminosities are at least one order 
of magnitude higher. We associate this discrepancy with their 
neglect of the nucleon-nucleon bremsstrahlung opacity and 
emission processes. Moreover, in the no-spin and co-rotating 
spin cases, the relationships between the fluxes of the var- 
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FIG. 10. — Colormaps of energy- and species-integrated specific neutrino 
energy deposition (whose volume-integrated values is refered to as Q(cc); 
see Fig. 17) and loss rates in the BNS merger model with no initial spin and 
shown here in units of 10 20 erg s _1 g -1 ). This snapshot corresponds to a time 
of 60 ms after the start of the simulation. The left section of the plot depicts 
the MGFLD result and the right shows the result of the S„ calculation. Note 
the distinctively enlarged polar gain regions and greater specific gain of the 
S„ result compared to the MGFLD calculation. This is in part a consequence 
of the larger polar neutrino fluxes and overall greater flux asymmetry in the 
S„ model (see Fig. 9). 

ious neutrino flavors changes significantly with time, partly 
because they are differently affected by neutron star cooling 
and changes in neutron richness. While the decrease of the 
luminosity predicted by VULCAN/2D makes physical sense, 
it contrasts with the findings of Setiawan et al. (2004, 2006) 
in their study of neutrino emission from torus-disks around 
~4M black holes. In their models, they incorporate phys- 
ical viscosity in the disk, together with the associated heat 
release, which keeps the gas temperature high. In our models, 
neutrino emission is a huge energy sink, not compensated by 
a-disk viscosity heating. 

The MGFLD treatment we use permits the calculation of 
the energy-dependent neutrino spectrum and its variation with 
angle. For the BNS merger models with initially no spins 
(left), co-rotating spins (middle), and counter-rotating spins 
(right), we show such spectra along the equatorial (black) and 
polar (red) directions in Fig. 8, using a reference radius of 
200 km. We can also explicitly compute average neutrino en- 
ergies, here defined as 



(4) 



J d£ v e 2 u F v (£ v ,R) 
J de v F v (e v ,R) 



(1) 



At 50 ms after the start of the VULCAN/2D simulations, 
such average neutrino energies in the no-spin BNS model and 
along the pole are 13.4 (v e ), 16.7 (i> e ), and 25.2 MeV ("V'), 
while they are 10.7, 15.6, and 2 1.2 MeV along the equatorial 
direction (in the same order). In the co-rotating model and 
along the polar direction, we have 14.0 (v e ), 17.0 (v e ), and 
22.0 MeV ("ffj,"), and along the equatorial direction, we have 
in the same order 9.9, 15.2, and 17.5 MeV. In the counter- 
rotating model and along the polar direction, we have 12.9 
(v e ), 17.1 (v e ), and 27.5 MeV ('V„"), and along the equatorial 
direction, we have in the same order 1 1 .5, 16.7, and 25. 1 MeV. 
The oblateness of the compact massive BNS merger, with 
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its lower temperature and its more gradual density fall-off 
along the equator, systematically places the neutrinospheres 
in lower temperature regions, translating into softer neutrino 
spectra. 

The aspherical density/temperature distributions naturally 
lead to a strong anisotropy of the radiation field, whose latitu- 
dinal dependence is shown in Fig. 9 for the three models (we 
use the same left to right ordering), and at a time of 60 ms after 
the start of the VULCAN/2D simulations. Note the stronger 
neutrino flux along the polar direction (see also Rosswog & 
Liebendorfer 2003, their Fig. 12), resulting from the larger 
radiating surface seen from higher latitudes, as well as the 
systematically larger matter temperatures at the decoupling 
region along the local vertical (also yielding a harder neutrino 
spectrum). In Fig. 9, we include in red the corresponding neu- 
trino fluxes computed with the 16 ^-angle S„ scheme, which 
shows comparable fluxes along the equator, but considerably 
larger ones at higher latitudes. We also show in Fig. 10 the 
energy- and species-integrated specific neutrino energy depo- 
sition and loss rates in the no-spin BNS model at 60 ms after 
the start of the simulation, for both the MGFLD and the S„ cal- 
culations. Notice how much more emphasized the anisotropy 
is with the multi-angle, S„, solver, and how enhanced is the 
magnitude of the deposition along the polar direction. This is 
a typical property seen for fast-rotating PNSs simulated with 
such a multi-angle solver (Ott et al. 2008). 

In Fig. 11, we show the energy-dependent neutrinosphere 
radii for the v e {left), D e {middle), and "v^" {right) neutri- 
nos, and for the BNS merger models with initially no spins 
{top row), co-rotating spins {middle row), and counter-rotating 
spins {bottom row). These decoupling radii manifest a strong 
variation with latitude, but also with energy due to the approx- 
imate e 2 dependence of the material opacity to neutrinos. The 
diffusion of neutrinos out of the opaque core, which leads to 
global cooling of the neutron star, is also energy dependent. 
The location-dependent diffusion timescale f* ff (r,z) is given 
by (Mihalas & Mihalas 1984; Ruffert et al. 1996): 

t^(r,z)^3T Vi AR/c, (2) 

where c is the speed of light, AR is the local density scale 
height, T Vi is the optical depth (species and energy dependent), 
and v\ is the neutrino species (i.e., v e , v e , or "i^"). For the v e 
neutrinos, we find an optical depth in the core that varies from 
a few x 10 2 for lOMeV neutrinos to a few x 10 4 for 100 Me V 
neutrinos. These optical depths translate into diffusion times 
that are on the order of 10 ms to 1 s. The quasi-Keplerian disk 
that extends from 20-30 to 100 km is moderately optically- 
thick to neutrinos at the peak of the energy distribution (i.e., 
at 10-20 Me V). Heat can thus leak out in the vertical direction 
over a typical diffusion time of a few tens of milliseconds (see 
Fig. 12). This is comparable to the accretion time scale of the 
disk itself (Setiawan et al. 2004, 2006), which suggests that 
this neutrino energy is available to power relativistic ejecta, 
and is not advected inwards with the gas, as proposed by Di 
Matteo et al. (2002). 

Compared to PNSs formed in the "standard" core-collapse 
of a massive star, SMNSs formed through coalescence have a 
number of contrasting properties: 1) most of the material is at 
nuclear density, unlike in PNSs where newly accreted mate- 
rial deposited at the PNS surface radiates its thermal energy as 
it contracts and cools; 2) most of the material is neutron rich, 
thereby diminishing the radiation of electron-type neutrinos 
in favor of antielectron-type neutrinos. In PNSs, electron- 



type neutrinos are an important cooling agent, but have large 
opacities. In SMNSs, the material is strongly deleptonized 
and electron-type neutrinos are non-degenerate. They suffer 
relatively lower opacity and can therefore diffuse out more 
easily. The opacity of other neutrino species is even lower 
and therefore more prone to energy leakage through radia- 
tion; 3) energy gain through accretion is modest since there 
is merely a few 0.1 M Q of material outside of the SMNS. In 
our three merger models, the internal energy budget of such 
shear-heated SMNSs is >10 53 erg, half of which is thermal 
energy. With a total neutrino luminosity of ^5 x 10 52 erg s -1 , 
the cooling time scale for these SMNS should be on the order 
of a second. This is significantly, but not dramatically, shorter 
than the 10-30 s cooling time scale of hot PNSs, i.e the time it 
takes these objects to radiate a total of ~3 x 10 53 erg of internal 
energy. 

A fraction of the diffusing core energy is absorbed through 
charge-current reactions (with associated volume-integrated 
energy £(cc)) and turned into internal energy in a "gain" layer 
at the surface of the SMNS, primarily on the pole-facing side 
and at high-latitudes (see Fig. 10). This neutrino energy de- 
position is the origin of a thermally-driven wind whose prop- 
erties we now describe. 

3.2. Neutrino-driven wind 

A few milliseconds after the start of our simulations and on- 
set of the burst of neutrinos from the BNS merger, a neutrino- 
driven wind develops. It relaxes into a quasi-steady state by 
the end of the simulations at 100 ms. Note that in this quasi 
steady-state configuration, and as shown in Fig. 10, the energy 
gain/losses due to neutrino in the polar direction is mostly net 
heating, with no presence of a sizable net cooling region. In 
Fig. 13, we show this evolution for the no-spin BNS merger 
model for four selected times: 10, 30, 60, and 100 ms. One 
can see the initial transient feature, which advects out and 
eventually leaves the grid. This transient phase is not caused 
by neutrino energy deposition, but is rather due to infall at 
small radii and in the polar regions of the ambient medium 
we placed around the merged object, as well as due to the 
sound waves generated by the core during the initial quakes. 
One can then see the strengthening neutrino-driven wind de- 
veloping in the polar funnel. The side lobes of the SMNS 
confine the ejecta at radii <500 km to a small angle of ^20° 
about the poles, but this opening angle grows at largers dis- 
tances to ^90°. The confinement at small radii also leads to 
the entrainment of material from the side lobes, overloading 
the wind and making it choke. The radial density /velocity 
profile of this wind flow is thus kinked, with variations in 
velocity that can be as large as the average asymptotic ve- 
locity value of ^30000 km s" 1 (Fig. 4). As shown in Fig. 13, 
the mass loss associated with the neutrino-driven wind is on 
the order of a few 10~ 3 Mq s" 1 str" 1 at a few tens of millisec- 
onds, but decreases to 10 -IO^Mq s" 1 str -1 at 100ms. The 
wind is also weaker along the poles than along the 70°-80° 
latitudes. The associated angle-integrated mass loss summed 
over 100 ms approaches 10~ 4 M Q and is made up in part of 
high F e material (~0.5) along the pole, but mostly of low Y e 
material (~0. 1-0.2) along mid-latitudes. Thus, "r-process ma- 
terial" will feed the interstellar medium through this neutrino- 
driven wind. This latitudinal variation of the electron fraction 
at large distances is controlled by the relative strength of the 
angle- and time-dependent v e and v e neutrino luminosities, 
the expansion timescale of the "wind" parcels, and the neutron 
richness at their launching site (Ott et al. 2007a). Along the 
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FIG . 11. — Using a colormap of the density distribution (shown in log scale) at 60 ms after the start of the VULCAN/2D simulations for the BNS merger models 
with initially no spins (top row), co-rotating spins (middle row), and counter-rotating spins (bottom row), we overplot in each case the contours corresponding to 
the energy-dependent and latitude-dependent neutrinosphere radii at neutrino energies e„ of 2.50, 6.87, 12.02, 21.01, 36.74, 64.25 MeV, for the u e (left column), 
P e (middle column), and "i/n," (right column) neutrinos. Corresponding radii (along a given latitude) increase monotonically with energy (matter opacity to 
neutrinos scales as ~ ej). For these calculations, the optical depth is integrated inwards along a fixed latitude starting at the maximum spherical radius of 
3000 km. 



pole and at the SMNS surface, the low-density, high neutron- 
richness, and relatively stronger v e neutrino luminosities at 
late times in the no-spin and co-rorotating spin models lead 
to a high asymptotic T e value (high proton-richness). De- 
spite the relatively high resolution employed in our simula- 
tions (^300 m in the radial direction at ^20 km), higher reso- 
lution would be needed to resolve this region accurately. Al- 
though we expect this trend would hold at higher resolution, 
it would likely yield lower asymptotic values of the electron 
fraction along the pole. 

Along the equatorial direction, low Y e material (~0. 1-0.2) 
migrates outward, but its velocity is below the local escape 
speed and it is unclear how much will eventually escape to 
infinity. 6 This is further illustrated in Fig. 14 where we show 

6 Note that in the counter-rotating model, an axis problem in the form of 



the angular variation of the density and velocity at 2900 km in 
the no-spin BNS model, at 121 ms. The BNS merger is thus 
cloaked along the poles by material with a density in excess of 
10 gem"" 3 , while along lower latitudes even denser material 
from the side lobes obstructs the view from the center of the 
SMNS. Importantly, wind material will feed the polar regions 
for as long as the merger remnant remains gravitationally sta- 
ble. Being so heavily baryon-loaded, the outflow can in no 
way be accelerated to relativistic speeds with high-Lorentz 
factors. In this context, the powering of a short-hard GRBs is 
impossible before black hole formation. 

a low-density, high-velocity, narrow region starts at the onset of the neutrino- 
driven wind and persists throughout the simulation. This spurious feature is, 
however, localized and therefore does not influence the global properties of 
the simulation. 
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FIG. 12. — Variation of the log of the optical depth along the z direction for 
the Ve neutrino at 12.02 MeV in the BNS merger model with no initial spins, 
plotted for a selection of cylindrical radii, 5, 15, 25, 50, and 100 km (see 
legend for the corresponding linestyles). The time corresponds to 60 ms after 
the start of the VULCAN/2D simulations. Notice that the disk transitions 
from optically-thick to optically-thin conditions at ~ 100 km, and that optical 
depth increases to a few hundreds only inside the SMNS. For most of the 
disk, the diffusion time (eq. 2) for 1 2.02 MeV u e neutrinos is on the order of 
a few milliseconds. The disk is therefore only moderately optically-thick for 
neutrinos with energies at the peak of the spectral energy distribution. 

As seen in Fig. 5, the rotational profile in the inner 100 km 
is strongly differential in the no-spin and counter-rotating spin 
cases, while it is quasi-uniform in the co-rotating spin case 
(see also Rosswog & Davies 2002, their Fig. 17). The good 
conservation of specific angular momentum in VULCAN/2D 
is in part responsible for the preservation of the initial rotation 
profile throughout the simulation. In reality, such a differen- 
tial rotation should not survive. Our 2D axisymmetric setup 
inhibits the development of tri-axial instabilities that arise at 
modest and large ratios of rotational and gravitational ener- 
gies (Rampp et al. 1998; Centrella et al. 2001; Saijo et al. 
2003; Ott et al. 2005, 2007b), i.e., under the conditions that 
prevail here. Moreover, our good, but not excellent, spatial 
resolution prevents the modeling of the magneto-rotational in- 
stability, whose effect is to efficiently redistribute angular mo- 
mentum (Balbus & Hawley 1991b; Pessah et al. 2006, 2007), 
and dissipate energy, and, in the present context, lead to mass 
accretion onto the SMNS. The magneto-rotational instabil- 
ity, operating on an rotational timescale, could lead to solid- 
body rotation within a few milliseconds in regions around 
~10km, and within a few tens of milliseconds in regions 
around ^-4 00 km. Note that this is a more relevant timescale 
than the typical ~100s that characterizes the angular mo- 
mentum loss through magnetic dipole radiation (Rosswog & 
Davies 2002). Importantly, in the three BNS merger models 
we study, we find that the free energy of rotation (the energy 
difference between the given differentially-rotating object and 
that of the equivalent solid-body rotating object with the same 
cumulative angular momentum) is very large. Modulo the 
differences of a factor of a few between models, it reaches 
~5xl0 51 erg inside the supermassive neutron star (regions 
with densities greater than 10 14 g cm" 3 ), but is on the order of 
2xl0 52 erg in the torus disk, regions with densities between 
10 11 and 10 14 gcm" 3 . Similar conditions in the core-collapse 
context yield powerful, magnetically- (and thermally-) driven 
explosions (LeBlanc & Wilson 1970; Bisnovatyi-Kogan et al. 
1976; Akiyama et al. 2003; Ardeljan et al. 2005; Moiseenko 



et al. 2006; Obergaulinger et al. 2006; Burrows et al. 2007a; 
Dessart et al. 2007). Rotation dramatically enhances the rate 
of mass ejection by increasing the density rather than the ve- 
locity of the flow, even possibly halting accretion and inhibit- 
ing the formation of a black hole (Dessart et al. 2008). In the 
present context, the magneto-rotational effects, which we do 
not include here, would considerably enhance the mass flux of 
the neutrino-driven wind. Importantly, the loss of differential 
rotational energy needed to facilitate the gravitational insta- 
bility is at the same time delaying it through the enhanced 
mass loss it induces. Work is needed to understand the sys- 
tematics of this interplay, and how much rotational energy the 
back hole is eventually endowed with. 

Oechslin et al. (2007), using a conformally-flat approxima- 
tion to GR and an SPH code, find that BNS mergers of the 
type discussed here and modeled with the Shen EOS, avoid 
the general-relativistic gravitational instability for many tens 
of milliseconds after the neutron stars first come into contact. 
Baumgarte et al. (2000), and more recently Morrison et al. 
(2004), Duez et al. (2004, 2006), and Shibata et al. (2006), 
using GR (and for some using a polytropic EOS), find that 
imposing even modest levels of differential rotation yields a 
significant increase by up to 50% in the maximum mass that 
can be supported stably, in particular pushing this value be- 
yond that of the merger remnant mass after coalescence. Sur- 
prisingly, Baiotti et al. (2008), using a full GR treatment but 
with a simplifed (and soft) EOS, find prompt black hole for- 
mation in such high mass progenitors. Despite this lack of 
consensus, the existence of neutron stars with a gravitational 
mass around 2 M Q favors a high incompressibility of nuclear 
matter, such as is in the Shen EOS, and suggests that SMNSs 
formed through BNS merger events may survive for tens of 
milliseconds before experiencing the general-relativistic grav- 
itational instability. In particular, the presence of a significant 
amount of material (a few tenths of a solar mass) located on 
wide orbits in a Keplerian disk, reduces the amount of mass 
that resides initially in the core, i.e., prior to outward transport 
of angular momentum. 

4. v-v ANNIHILATION 

Besides being prime candidates for gravitational wave 
emission, binary neutron star mergers may lead to short-hard 
GRBs (for an overview, see, e.g., Piran 2005; Nakar 2007). In 
this context, the high energy radiation arises from non-thermal 
radiation associated with relativistic ejecta. The annihilation 
of neutrino pairs, by the process 

Vj + v, ^ e + + e~ ; i e {e, fx, r} , (3) 
represents one possible powering source for high-Lorentz 
factor baryon-free ejecta beamed into a small solid angle 
about the rotation axis of the BNS merger (Ruffert et al. 1997; 
Ruffert & Janka 1999; Asano & Fukuyama 2000, 2001 ; Miller 
et al. 2003; Kneller et al. 2006; Birkl et al. 2007; Rosswog 
et al. 2003; Rosswog 2005, and references therein). 

Here, we revisit this proposition by studying quantitatively 
the energetics of neutrino-antineutrino annihilation in the 
three BNS mergers we simulated with VULCAN/2D. First, 
in § 4.1, we apply the approach used in previous work and 
based on the leakage scheme of Ruffert et al. (1996, 1997). 
This approach was also used by Rosswog & Liebendorfer 
(2003), but for merger events whose dynamical evolution was 
computed independently. 7 Note that similar VjVj annihilation 

7 SPH versus grid-based code; different EOS; different resolution. 
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FIG. 13. — Colormaps of the log of the mass loss rate per steradian (d 2 M/dtdCl, in units of Mq s ' str ) for the no-spin BNS merger model at 10ms (top 
left), 30 ms (top right), 60 ms (bottom left), and 100 ms (bottom right) after the start of the VULCAN/2D simulation, and depicting the mass loss associated with 
the initial transient, followed by the neutrino-driven wind. The displayed region covers 2000x2000 km 2 . Regions that are infalling or denser than 10 10 gcm~ 3 
are shown in red, and velocity vectors, overplotted in black, have a length saturated at 7% of the width of the display for a magnitude of 30,000 kms , Notice 
the concomitant mass loss from the poles down to mid-latitudes (the wind) and the expansion of BNS merger material at near-equatorial latitudes. 



rate calculations, i.e., without the full momentum-space an- 
gular dependence, have been carried out in the core-collapse 
and PNS contexts, both with Newtonian gravity (Goodman 
et al. 1987; Cooperstein et al. 1986, 1987; Janka 1991) and in 
GR (Salmonson & Wilson 1999; Bhattacharyya et al. 2007). 
Then, in § 4.2 we present a new formalism based on moments 
of the neutrino specific intensity. The angle and energy de- 
pendence of the neutrino radiation field is obtained by post- 
processing individual MGFLD VULCAN/2D snapshots (see 
§ 3) with the multi-angle, S„, variant. In both cases, because 
the power associated with neutrino-antineutrino annihilation 
is sub-dominant compared to that associated with the charge- 
current reactions prior to black hole formation, it can be esti- 
mated only through a post-processing step. 

4.1. Formalism based on a leakage scheme and the approach 
ofRuffertetal. (1996) 



Estimation of the neutrino-antineutrino annihilation rate 
using the Ruffert et al. (1996, 1997) approach is done in 
two steps. First, using the processes described in §3.1, 
the instantaneous rate of neutrino emission Q(yD is com- 
puted for all grid cells. It is subsequently weighted by the 
direction-dependent factor that depends on the radiative dif- 
fusion timescale ?* ff and neutrino emission timescale t x ° ss rel- 
evant for that cell. Since we are primarily interested in the 
energy deposition in the polar regions, we consider only the 
cylindrical-z direction when estimating fj™. The effective 
emissivity Q eS (fi) then has the form 



1 + (f£ ss ) t*f 

where expressions for the various components are given ex- 
plicitly in Appendices A & B of Ruffert et al. (1996). 

In Fig. 15, we show the effective emissivity resulting from 
this leakage scheme and with the neutrino processes of Ruf- 



14 



TABLE 1 

Ujpj ANNIHILATION RATES USING THE LEAKAGE AND THE S„ SCHEMES 







No Spins 


Co-Rotatinj 


> Spins 


Counter-Rotating Spins 


Time 


J dVQ + (u e D e ) 




}dVQ + (u e D e ) 






J' dV(F(y e D e ) 




(ms) 




Bs-' 




Bs 1 






Bs" 1 


Leakage Scheme 


10 


1.56(-1) 


3.28(-5) 


2.05(-l) 




1.01(-4) 


2.18(-1) 


1.40(-4) 


20 


9.28(-2) 


5.99(-6) 


4.79(-2) 




2.38(-6) 


1.05(-1) 


1.01 (-5) 


30 


3.80(-2) 


1.15(-6) 


3.70(-2) 




1.43(-6) 


5.95(-2) 


3.35(-6) 


40 


3.11(-2) 


8.54(-7) 


4.67(-2) 




1.08(-6) 


3.65(-2) 


1.36(-6) 


50 


2.82(-2) 


6.23(-7) 


5.67(-2) 




1.35(-6) 


2.16(-2) 


5.81(-7) 


60 


1.78(-2) 


2.57(-7) 


4.43(-2) 




8.74(-7) 


1.50(-2) 


3.00(-7) 


70 


1.32(-2) 


1.63(-7) 


2.69(-2) 




3.78(-7) 


1.18(-2) 


1.79(-7) 


80 


1.22(-2) 


1.67(-7) 


3.55(-2) 




8.16(-7) 


1.01 (-2) 


1.21(-7) 


90 


1.47(-2) 


2.94(-7) 


2.96(-2) 




6.64(-7) 


7.64(-3) 


7.64(-8) 


100 


1.60(-2) 


3.73(-7) 


2.54(-2) 




5.41(-7) 


6.48(-3) 


5.47(-7) 


S n Scheme 


10 


1.81(-1) 


1.44(-2) 


5.97(-l) 




9.19(-3) 


3.63(-l) 


4.76(-2) 


20 


6.82(-2) 


6.41(-3) 


9.22(-2) 




1.15(-3) 


8.92(-2) 


1.70(-2) 


30 


3.95(-2) 


3.96(-3) 


4.59(-2) 




4.59(-4) 


4.08(-2) 


1.01 (-2) 


40 


2.71(-2) 


2.58(-3) 


4.13(-2) 




2.82(-4) 


2.77(-2) 


7.83(-3) 


50 


2.18(-2) 


1.78(-3) 


4.19(-2) 




1.89(-4) 


1.92(-2) 


5.84(-3) 


60 


1.82(-2) 


1.39(-3) 


3.59(-2) 




1.34(-4) 


1.28(-2) 


4.27(-3) 


70 


1.47(-2) 


1.02(-3) 


2.31(-2) 




8.70(-5) 


1.04(-2) 


3.76(-3) 


80 


l.ll(-2) 


6.94(-4) 


2.30(-2) 




6.22(-5) 


8.49(-3) 


3.28(-3) 


90 


1.01(-2) 


5.25(-4) 


1.84(-2) 




4.28(-5) 


7.13(-3) 


2.83(-3) 


100 


9.67(-3) 


3.93(-4) 


1.59(-2) 




3.02(-5) 


6.21(-3) 


2.46(-3) 



NOTE. — Upper half: Listing of the volume -integrated rates for the i/ e V e and "f „ annihilation rates using the approach of Ruffert et al. (1996, 1997), as summarized in §4. 1. 
The time for each calculation refers to the time since the start of the VULCAN/2D simulation. Numbers in parenthesis correspond to powers of ten. The strong decrease of the energy 
deposition rates reflects the fading of the neutrino luminosity due to the cooling of the SMNS. A density cut above 10 1 1 g cm -3 is applied for both emission and deposition sites. Lower 
half: Same as upper half, but this time using the 8 $-angle S„ calculations, performed through a post-processing of the MGFLD radiation -hydrodynamics snapshots computed with 
VULCAN/2D. The same hydrodynamics background (p, 7", Y c distribution) is used for the leakage and the -schemes. 
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FIG. 14. — Variation of the radial velocity (solid) and density (dotted; Log 
scale) versus colatitude and at a spherical radius of 2900 km for the BNS 
merger model with no initial spin. The time plotted is 121 ms after the start 
of the MGFLD VULCAN/2D simulation. Note that the escape speed at this 
radius and for a ~2.8Mq central object is ~16000 kms~', so most of the 
mass near equatorial latitudes will remain trapped in the gravitational poten- 
tial of the SMNS. 



fert et al. (1996) for the v e {left column), v e (middle col- 
umn), and "i/„" {right column) neutrinos, for the BNS merger 
models with initially no spins {top row), co-rotating spins 
{middle row), and counter-rotating spins {bottom row). For 
each process, emission is conditioned by the competing el- 
ements of optical depth, density, temperature, and electron 
fraction. In practice, it peaks in regions that are dense, hot, 
but not too optically thick, i.e., at the surface of the SMNS 
(see density contours in Fig. 15, overplotted in black). Op- 
tical depths in excess of 100 for all neutrino energies make 
the inner region (the inner ~15 km, where densities have nu- 
clear values) a weak "effective" emitter, their contribution 
operating on a 0. 1—1 s timescale. The dominant emission 
associated with the neutrinos originates from a con- 



siderably higher-density region than that associated with the 
electron-type neutrinos, i.e. 10 12 -10 13 g cm" 3 compared to 
10 9 -10 n gem" 3 , with a corresponding "leakage" luminosity 
in all three merger models that is typically an order of mag- 
nitude smaller than predicted by VULCAN/2D (see Fig. 7). 
This low "Up" emissivity is likely caused by the neglect of 
nucleon-nucleon bremsstrahlung processes in the approach of 
Ruffert et al., which leads to a smaller decoupling radius for 
"Vu" neutrinos, and, therefore, an underestimate of the size 
of the radiating surface from which they emerge (shown in 
Fig. 11, right column). In practice, the decoupling radius is 
energy-dependent, but here we stress the systematic reduc- 
tion of the decoupling radius for all 'V M " neutrino energies 
because of the neglect of this extra opacity source. The im- 
portance of the bremsstrahlung process for "v^' emissivity 
and opacity has been emphasized by Thompson et al. (2000) 
for "hot" PNSs. The SMNS that results from the merger is 
considerably heated by shocks and shear, and is thus also in 
a configuration where such bremsstrahlung processes are im- 
portant and should be included. 

From the effective emissivity distribution computed for 
each neutrino flavor with the leakage scheme, Ruffert et al. 
(1997) then sum the contributions from all pairings between 
grid cells. For completeness, we briefly reproduce here the 
presentation of Ruffert et al. (1997). The integral to be com- 
puted for the energy-integrated z/,i/; (representing equivalently 
v e v e or "VfiVfx") annihilation rate at position r is 
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FIG. 15. — Colormap of the log of the "effective" emissivity Q eff (Ui) for u e (left), u e (middle), and "uu" (right) neutrinos at 60ms after the start of the 
simulation for the BNS merger configuration with no initial spins (top row), co-rotating spins (middle row), and counter-rotating spins (bottom row). We also 
overplot iso-density contours at every decade between 10 10 and 10 14 gem -3 , to allow a visual assessment of the effect of the adopted density cut (at 10" gem -3 ) 
on the computation of the annihilation rate. Note how severely this cut truncates the emission of "i/^" neutrinos. The calculation is based on the formalism 
developed and presented by Ruffert et al. (1996). 

as § 4.2). (To is the baseline weak interaction cross section, 
1.705 xl0~ 44 cm 2 , c is the speed of light, m e is the electron 
mass. <f> is the angle between the neutrino and antineutrino 
beams, entering the annihilation rate formulation through the 
term (l-cos<I>) (squared or not) which thus gives a stronger 
weighting to larger angle collisions. This is what makes the 
dumbell-like morphology of BNS mergers such a prime can- 
didate for v-,Vi annihilation over spherical configurations (see 
also § A). (e) Vj and (e) Pf are the Vi and Vi mean neutrino ener- 
gies, whose values we adopt for consistency from the simula- 
tions of Ruffert et al. (1997), i.e., 12, 20, and 27MeV for u e , 
D e , and "v^", respectively (our values are within a few times 
10% of these, so this has little impact on our discussion). Note 
also that the average neutrino energies for all three species are 
much larger than m e c 2 , and therefore make the second term 
in the above equation negligible (it is about a factor of 1000 



„+. 1 oo f(Ci+C 2 ) 1/P 

Q wvo = ^ ■ 



I 6£tI Vi <f dfi'ip, [(e) Vi + (e) 9l ] (l-cos$) 2 + 

J Air J At: 



+Cj,, Vi 9. (m e c 2 y 



I &ll vt I dfi% f fc\"' + ^ p ' (l-cos$) 1 . (5) 

J Air J Air (^i(e)Pi J 

I Vi is the neutrino specific intensity, ft and S7' are the solid 
angles subtended by the cells producing the neutrino and an- 
tineutrino radiation incident from all directions. C\, C2, and 
C3 are related to the weak coupling constants Ca and Cy and 
depend on the neutrino species (see Ruffert et al. 1997 as well 
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smaller than the first one, and also has a much weaker large- 
angle weighting). The total annihilation rate is then the sum 

G + (^ e )+G + ("ivV)- 

In practice, we turn the angle integrals into discretized sums 
through the transformation 



dfLL, 



k 



k -lv,k , 



(6) 



where Af2* is the solid angle subtended by the cell k as seen 
from the location 7. Ruffert et al. (1997) applied their method 
to 3D Cartesian simulations, while our simulations are in 2D 
and have both axial symmetry, and mirror symmetry about the 
equatorial plane. Making use of these symmetry properties, 
we remap our VULCAN/2D simulation from a 2D (cylindri- 
cal) meridional slice onto a 3D spherical volume which covers 
the space with a uniform grid of 40 zones for the 2tt azimuthal 
direction and 20 zones for the tt polar direction, while the re- 
duced radial grid has a contant spacing in the log and uses 20 
zones between 12 and 120 km. We compute the annihilation 
rate at locations in a 2D meridional slice of this new volume, 
with 16 uniformly-spaced angles from the pole to the equator, 
and 40 zones with a constant spacing in the log between 15 
and 120 km. 8 To estimate the annihilation rate at r, one needs 
to estimate the flux received from all cells r^. Ruffert et al. 
(1997) assume that neutrino radiation is isotropic in the half 
space around the outward direction given by the density gra- 
dient n p at k. Moreover, they approximate each emitting cell 
volume as a sphere, whose radius at (r k , 9 k , 4>k) is 



Dv = 



4/r 



(7) 



where fj, k = cosf^. With our location dependent cell volume, 
we obtain 



k 



Ml k I v i 



-£^#^V,), 



34 



(8) 



where du = |r — r^|. When we sum over discrete cells, the 
cosine of the angle between k and k' is given by cos <&kk' = 
if- n)ir- r k i)/i\r- r k \\r- f%> |). 

Following Ruffert et al. (1997), we apply selection crite- 
ria to determine whether to include a contribution. Specifi- 
cally, emission and deposition sites must be in regions with 
a density lower than 10" gem" 3 . Interestingly, the leakage 
scheme predicts a very small "effective" emissivity from high- 
density regions, due to the very long diffusion times from 
those optically-thick regions. Relaxing the density cuts in 
the calculations leads only to a 50% enhancement in volume- 
integrated annihilation rate E(v e D e ), but an increase of a factor 
of ~20 for ECu^v^,"). As pointed out earlier, the location in 
very high-density regions of the 'V^" emitting cells means 
that most of the emitting volume is truncated by the adopted 
10 11 gem -3 density cut, while results converge when this cut 
is increased to densities of ~10 13 gem -3 . Even with the lat- 
ter, EC'VfiVn") is still three orders of magnitude smaller than 
£(^ e P e ), likely a result of the neglect of bremsstrahlung, and 
the unfavorable (1 — cos$^/) 2 for this compact emitting con- 
figuration. 

By contrast with the rather large uncertainty introduced 
through the somewhat arbitrary density cuts, the annihila- 
tion rate calculation is weakly affected by increasing the 

8 We also perform higher resolution, accuracy-check, calculations, but our 
results are already converged at this lower resolution. 



resolution. In the no-spin BNS model at 50 ms, chang- 
ing the number of radial-angle zones for the deposition 
sites {«r ec je P ,rcfedep} from (40,16) to (60,48) and for emit- 
ting sites {nr em it,«Wit,"Pemit} from (20,20,40) to (30,30,60) 
increases E(v e D e ) from 2.82xl0 49 to 3.14xl0 49 ergs _1 (a 
~10% increase) and increases E^v^v^") from 6.23 x 10 44 to 
7.59 x 10 44 erg s" 1 (a ~20% increase). The same test done for 
the no-spin BNS model at 60 ms after the start of the simula- 
tions yields a E(v e D e ) which is identical to within 1%, while 
Ei"VuVu") increases by ~10% in the higher resolution model. 
Hence, higher resolution does not change these values by a 
significant amount. 

In the left half of each panel in Fig. 16, we show the dis- 
tribution in a 2D meridional slice (it is in fact axisymmetric 
by construction) for the annihilation rate Q + (v e v e ) of electron- 
type neutrinos, for the BNS merger models with initially no 
spins (left), co-rotating spins (middle), and counter-rotating 
spins (right). Volume-integral values E(v{Vi) as a function of 
time for all three models are shown in Fig. 17 (left panel; 
dashed line joined by star symbols) and given in Table 1 . The 
deposition rate is maximum near the poles, where the large- 
angle collisions occur, and close to the peak emission sites 
(shown in Fig. 15), i.e., near the SMNS surface, with peak 
values on the order of 10 29 ergcm~ 3 s" 1 for all three mod- 
els. In this formalism, Ei^v^v^) (middle panel of Fig. 17) 
is at least three orders of magnitude smaller than E(v e D e ). 
The peak volume-integrated energy deposition reaches a few 
10 50 ergs _1 , typically two orders of magnitude below that 
achieved by charge-current reactions, whose associated en- 
ergy deposition drives the neutrino-driven wind (§ 3; right 
panel of Fig. 17). 

4.2. Formalism based on multi-angle transport and moments 
of the specific intensity 

We now present a formalism for the computation of the i/,z/,- 
annihilation rate that directly exploits the neutrino-transport 
solution computed with VULCAN/2D in our BNS mergers, 
rather than using evaluations based on the leakage scheme 
(Ruffert et al. 1996, 1997). 

Our approach is to use the angle-dependent neutrino spe- 
cific intensity calculated for snapshots at 10 ms intervals for 
these BNS merger models using the multi-angle, S„, scheme 
(Livne et al. 2004; Ott et al. 2008; § 2.2). In the Appendix, 
and for completeness, we present such annihilation rate cal- 
culations, but in the context of the post-bounce phase of a 
20 M Q progenitor (Ott et al. 2008) and of the accretion in- 
duced collapse (AIC) of a massive and fast rotating white 
dwarf (Dessart et al. 2006b). The formalism presented here 
applies to all cases equivalently. 

Following Burrows et al. (2006) and Janka (1991), the ex- 
pression for the local energy deposition rate Q + (vjh>i) at a po- 
sition Fdue to annihilation of v(Di pairs into electron-positron 
pairs is given to leading order by 

1 cra(C\ + Cy) Vii>i 



Q + im)= 



6 c(m e c 2 ) 2 



pOO poo 

■ / de v / ds' v (s v + e'j x 
Jo Jo 



I dQf dO'V^l-cos*) 2 , (9) 

J Ait J Ait 

where l Vl = I Vi (e v , r, n, t) is the z^-neutrino specific intensity at 
energy e v , location r, along the direction n, and at time t. The 
primes denote the anti-particle. The angle $ is the angle be- 
tween the directions n and ff of the neutrino and antineutrino, 
i.e., cos$ = n-fC. 
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FIG. 16. — Colormaps of the log of the energy deposition Q + (u e u s: ) by electron neutrino pair annihilation in the BNS merger models with initially no spin (left), 
co-rotating spins (middle), and counter-rotating spins (right), and at 60 ms after the start of the VULCAN/2D simulation. In each case, the evaluation is done with 
two different methods. In the left half of each panel, we show the results using the leakage scheme (Ruffert et al. 1996, 1997; see also § 4.1). In the right half, 
we show the corresponding results using the new formalism presented in § 4.2 using 16 ^-angles. We set the minimum of the colorbar to white for all regions 
with densities greater than 10' 1 gcirr 3 . The maximum of the color bar is set at 10 29 ergenr 3 s to improve the visibility. Maximum values differ in fact by 
large factors between the two methods (in contrast with the good agreement in the cumulative values). For the model with no initial spins, we obtain maxima at 
l.Ox 10 30 (S„), compared with 6.5 X 10 28 erg cirT 3 s _I (leakage), and in the same order, we have 5.2 X 10 29 and 1.4 X 10 29 erg cm~ 3 s for the co-rotating model, 
and 1.4xl0 30 and 2.1 xlO 29 erg cm s for the counter-rotating model. We do not show the distribution for the "v^v^" annihilation, which is qualitatively 
similar to that of v t v s , being merely weaker everywhere by about an order of magnitude. With the leakage scheme, the contrast in annihilation rate between 
these two neutrino types is even greater, primarily because 'V M " neutrino emission occurs at densities in excess of 10 12 gcnT 3 , thus, beyond our density cut. 
In the snapshots shown here and with the leakage scheme, we find integral net energy deposition rates by neutrino pair annihilation E(v\v\) (summed over all 
neutrino species) of 1.78x 10 49 , 4.43x 10 49 , and 1.50x 10 49 erg s _1 for the BNS merger models with initially no spins, co-rotating spins, and counter-rotating 
spins, respectively. In the same order, but with the S„ scheme (and using 16 i?-angles), we find 2.13x 10 49 , 3.54 X 10 49 , and 1.44x 10 49 erg s _I . Note that these 
numbers are at least one to two orders of magnitude smaller than the corresponding rates due to charged-current neutrino absorption (see left and right panels of 
Fig. 17). 



The formulation we present applies equally to all neutrino 
species, as long as the values of the weak coupling constants 
Ca and Cy are appropriately set. We have Cy = 1 /2 + 2 sin 2 <9 W 
for the electron types, Cy = -1 /2 +2 sin 2 f?w for the and v T 
types, and Ca = ±1 /2, with sin 2 #w = 0.23. The other variables 
have their usual meanings. Note that this formulation neglects 
the phase-space blocking by the final-state electron-positron 
pair which is relevant only at high densities and negligible 
in the semi-optically thin regime where pair annihilation may 
contribute to the net heating. 

Now, setting the flux factor h = H/J and the Eddington- 
tensor factor k = K/J and by expanding the cos$ = ft ■ ft 
term with an appropriate choice of the radiation unit vector 
(Hubeny & Burrows 2007), we obtain for the general case in 
three dimensions 



io^jCl + Cy)^ 
3c(m e c 2 ) 2 

J Vl j'(l-2h Ul -h' +Tr[k„,:k' ]) , 



pOO pOO 

/ de v \ de' v (e v + e' v ) x 
Jo Jo 



(10) 



where the term Tr[k„, : k p .] is the trace of the matrix product 
of the two J u -normalized Eddington tensors. 

Assuming the special choice of cylindrical coordinates in 
axisymmetry and neglecting the velocity dependence of the 
radiation field, we obtain 



Za^iCl + Cy)^, 

3c(m e c 2 ) 2 
J J 1 x(l-2h r h' r -2h z h' z 

+KX rr + k, z k^ + k^k^ + 2k r ,k^) 



pOC pOO 

/ de v / de' v (e v + e' v ) x 
Jo Jo 



(11) 



where we have suppressed the v-, and v, subscripts. Note that 
k rz is the only nonzero off-diagonal term in these coordinates. 
Using the trace condition on the Eddington tensor, k rr + W zz + 



= 1, the above can be recast into 



Q + (Wi) = 



pOO pOC 

/ de v / de' y {£ u + e' u ) x 
Jo Jo 



8o~07T (Cfl+Cy)^,. 

3c(m e c 2 ) 2 
J J' x ( 1 - 2h r h' r - 2h z h' z + W rr V! rr + K- ; k', 

+( 1 - K, - k„)( 1 - Kr ~ Kz) + 2 krX z ) ■ 



(12) 



) X 

(13) 



For completeness and comparison with Janka (1991), we 
transform to spherical coordinates and reduce to spherically 
symmetric problems, which then yields 

, _ 8(7o7r 2 i o f°° f°° 

& ! (w)=T7 2^( C A +C v)m / d£ » / de' v {e v + e' v 

ic(m e c ) Jo Jo 

// x (\-2hh'+kk' + ^{\-k){\-k')) , 

with scalar h and k. Equation (13) corresponds to the combi- 
nation of eqs. (1) and (6) of Janka (1991). 

In practice, our approach is to start from snapshots of the 
(converged) MGFLD radiation-hydrodynamics simulations of 
the three BNS mergers presented in § 3. Keeping the hydro- 
dynamical variables frozen, the radiation variables are relaxed 
using the S„ algorithm and with 8 ^-angles. As an accuracy 
check, we also compare the S„ results based on simulations 
performed with 16-i9 angles. When we do this, the converged 
radiation field for the 8 $-angle solution is remapped into 
the 16 i?-angle run (corresponding to a total of 144 angles), 
which is then relaxed. These simulations are quite costly, and 
take about 4 days with 48 processors. However, by contrast 
with the approach of Ruffert et al., the annihilation rates are 
then simply obtained from direct (local) integration of the mo- 
ments of the neutrino specific intensity. 

In the right half of each panel in Fig. 16, we show the dis- 
tribution in a 2D meridional slice for the annihilation rate 
Q + (v e v e ) of electron-type neutrinos, for the BNS merger mod- 
els with initially no spins (left), co-rotating spins (middle), 
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FIG. 17. — Volume-integrated energy deposition rate E(v c v c ) (top) and the 
EC'i/fiD^") (middle) neutrino-antineutrino annihilation processes, as well as 
that due to charge-current reactions (bottom) for the models initially with no 
spins (black), co-rotating spins (red), and counter-rotating spins (blue) and 
computed with the formalism presented in § 4.2. The dashed lines (symbols 
highlight the times computed) represent for each model the corresponding 
result based on the leakage scheme (§ 4.1), following the method of Ruf- 
fert et al. (1996, 1997). The computations are performed every 10 ms, from 
10 to 100 ms after the start of the VULCAN/2D simulations. Note the fast 
decrease of the energy deposition rate with time (oc f" 1 8 ; dotted line), follow- 
ing the decrease of all neutrino luminosities (see Fig. 7) and the contraction 
of the SMNS. Note that Setiawan et al. (2004, 2006) obtain a slightly weaker 
dependence (oc f" 15 ), but with a black-hole/torus-disk configuration and a 
treatment of the physical viscosity for energy dissipation through shear in the 
differentially-rotating disk. Note that the energy unit used is the Bethe, i.e. 
10 51 erg= 1 Bethe [IB]. 

and counter-rotating spins (right), but now computed from 
the results of the S n calculation and the formalism presented 
above. We provide in Fig. 17 volume-integral values as a 
function of time for all three models {top and middle pan- 
els; solid lines) and given in Table 1. By contrast with the 
leakage scheme results, the deposition is more evenly spread 
around the SMNS, is maximum at depth (in the region where 
the density cut applies), and reaches peak values a factor of 



a few larger. We find that it is not only the large-angle col- 
lisions that favor the deposition, but also the dependence on 
J 2 (and its l/R 4 radial dependence in optically-thin regions) 
which considerably weights the regions close to the radiat- 
ing SMNS, making the deposition large not only in the polar 
regions, but also all around the SMNS. After 40ms, and al- 
though the associated neutrino luminosities are quite compa- 
rable between the three models, the (somewhat unrealistic) 
co-rotating spin BNS model boasts the largest annihilation 
rate. This stems from the very extended high-density SMNS 
configuration, leading to larger neutrinosphere radii and ex- 
tended regions favoring large-angle collisions. This constrasts 
with the fact that it is the least hot of all three configurations, 
and is a somewhat weaker emitter, suggesting that it is not 
just the neutrino luminosity, but also the spatial configuration 
of the SMNS that sets the magnitude of the annihilation rate. 
In the new formalism based on the S„ simulations, E^'v^v^") 
{middle panel of Fig. 17) is now only one order of magnitude 
smaller than E{v e i> e ), even for the same adopted density cut. 
By contrast with the leakage scheme, and to a large extent be- 
cause of the treatment of bremsstrahlung processes, the emis- 
sion from "Vft" neutrinos is associated with more exterior re- 
gions of the SMNS, the ~25 MeV "v^' neutrinos decoupling 
at ~ 1 8 and ~ 1 00 km along the polar and equatorial directions, 
respectively, and thus in regions with densities on the order of, 
or lower than, 10 11 g cm" 3 (Fig. 11). The emission is, thus, not 
significantly truncated by the adopted density cut. Moreover, 
the subtended angle of the representative 'V M " neutrinosphere 
is correspondingly much bigger, favoring larger-angle colli- 
sions. 

This disagreement should not overshadow the good match 
between the leakage- and the ^-scheme predictions for the 
i/ e D e annihilation rate, which are typically within a few tens 
of percent, only sometimes in disagreement by a factor of 
~3 (see, e.g., the no-spin BNS model at 10 ms). The S n - 
scheme results, thus, support the long-term decline of the an- 
nihilation rate predicted with the leakage scheme, but with 
a EC'i^iy,") value that is typically a factor of ten smaller 
than E{i/ e i> e ) at all times. The results given here using 8 
^-angles are already fairly converged. For the co-rotating 
spin BNS model at 60 ms after the start of the simulation, 
we obtain the following differences between the 16 and 
8 i?-angle S n calculations: E{i/ e i> e ) changes to 3.54 xlO 49 
from 3.59 x 10 49 erg s" 1 (down by 2%), E^v^v^) changes to 
1.73xl0 47 from 1.34x 10 47 erg s" 1 (up by 29%), and £(cc) 
changes to 1.40x 10 51 from 1.64x 10 51 ergs" 1 (down by 15%) 

The long time coverage of our simulations, up to > 100 ms, 
shows that the peak value, which is also coincident with the 
peak neutrino luminosity at 5-10 ms, is not sustained. The 
conditions at the peak of the annihilation rate do not cor- 
respond to a steady-state, but instead herald the steady de- 
crease that follows as the BNS merger radiates, contracts, and 
cools. In all simulations, E{v e v s ) is down by a factor 10-100 
at 100 ms compared to its peak value, while for E^'v^v^") the 
decrease is by 2-3 orders of magnitude (this component is in 
any case sub-dominant). The J 2 dependence, combined with 
the strong fading of all neutrino emissivities, is at the origin 
of the steady decrease of the annihilation rate over the time 
span considered here. In their black-hole/torus-disk config- 
uration with a-disk viscosity (which causes significant shear 
heating not accounted for here), Setiawan et al. (2004, 2006) 
observe that the annihilation rate decreases as f" 3 ' 2 , and, thus, 
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FIG. 18. — Angular factor W in the VjVj annihilation rate (see eqs. 11 and 
14) shown for the v e v e process at e„ = 12.02 MeV for the BNS model with 
no initial spin. The snapshot corresponds to a time of 60 ms after the start of 
the VULCAN/2D simulation. We superpose v e flux vectors at the same &v 
and density contours corresponding to 10 9 , 10 10 , and 10" gem -3 , indicating 
the local density gradient with perpendicular tick marks. 

only slightly more gradually than our prediction. In our work, 
neutrino energy deposition by charge-current reactions is the 
dominant means to counteract the global cooling effect of 
neutrino emission. 

We now use the results for the S n scheme to discuss the 
assumption of Ruffert et al. that neutrino emission is quasi- 
isotropic. More precisely, they propose that the effective 
emission for every cell is isotropic in the half space around 
the outward direction given by the density gradient n p . By 
contrast, in the S„ scheme, the angular distribution of the neu- 
trino radiation field is solved for, and we can determine how 
isotropic this emission is. For that purpose, we study the spa- 
tial variation of the angular term in eq. (11), 

W = (l-2h,.h' r -2h z h' z + 

and plot it in Fig. 18 for the v e —v e annihilation rate at a rep- 
resentative neutrino energy of 12.02 MeV. Overplotting the 
(representative) 12.02 MeV v e flux vectors, together with iso- 
density contours, one sees that the flux is oriented everywhere 
predominantly in the radial direction, and peaks along the po- 
lar direction. The outward direction given by the local density 
gradient (given by the perpendicular to iso-density contours, 
shown here as short white tick marks), is colinear to the flux 
vector along near-polar latitudes, while at mid-latitudes, they 
are perpendicular to each other. Hence, the assumption by 
Ruffert et al. of isotropic neutrino emission in the half space 
around the outward direction given by the local density gra- 
dient is not accurate in this instance. Note that it will likely 
hold more suitably once the black hole forms and only the 
torus disk radiates. 

In Fig. 18, W has a maximum of ~1.33 in the optically- 
thick regions of the extended (equatorial) disk-like structure. 
There, the flux-factor terms and k' w + 2k rz k' r , are nearly zero, 
while k rr k' rr + k zz k' zz + k vv k' vv = 1/3. In polar regions, the tran- 
sition to free streaming happens at small radii and W de- 
creases from ~1 .33 at <20km to ^0.3 at 40 km. Outside 
~60-100km, W w 0, since h zz h' zz w k zz k' zz w 1 and all other 



terms are small. Considering now the fact that the distribu- 
tion of the mean intensity, J, irrespective of energy group 
and species, is oblate inside the rapidly spinning supermas- 
sive neutron star and transitions to a prolate shape outside (see 
Fig. 10), we can explain the spatial distribution of the pair an- 
nihilation rate (as shown in Fig. 16). The rate peaks inside the 
oblate SMNS, since J J 1 and W are largest there. Along the po- 
lar axis, W decreases rapidly with radius, while //' becomes 
prolate and remains large out to large radii, compensating in 
part for the small W. In equatorial regions, on the other hand, 
J J' drops off more rapidly, but W remains near 1.33 out to a 
radius of ^150 km, and only slowly transitions to zero as the 
equatorial radiation field gradually becomes forward-peaked. 
This systematic behavior sustains a significant V(Di annihila- 
tion rate for equatorial radii of up to ~ 1 50 km. Note, however, 
that net energy deposition by pair annihilation does not occur 
in the equatorial charged-current loss regions. 

The findings described in the above paragraph, indicate 
that the dominant contribution to Q + (vii>i) is not coming from 
large-angle collisions of neutrinos emitted from the disk, but 
rather is due to collisions at all angles and close to the high- 
density high-temperature SMNS surface. This is in distinction 
to the annihilation rate computed using the leakage scheme 
and the assumptions concerning the morphology of the neu- 
trino radiation field made in Ruffert et al. (1997). However, 
once the SMNS transitions to a black hole, the inner contribu- 
tion will vanish and the radiation will come exclusively from 
the cooling hot torus. We will address in a future paper what 
annihilation rates we obtain with this BNS merger configura- 
tion. 

In the future, while retaining the merits of the S„ method, 
we will need to improve the computation of the annihilation 
rates by accounting for GR effects. The compact and massive 
configuration of BNS mergers, with GM/Rc 2 on the order of 
20% at 20 km, suggests that the numbers we present could be 
modified with this improvement, although Birkl et al. (2007) 
suggest the magnitude of the effect is at most a few tens of 
percent. 

5. CONCLUSION 

We have presented multi-group flux-limited-diffusion radi- 
ation hydrodynamics simulations of binary neutron star merg- 
ers, starting from azimuthal-averaged 2D slices based on the 
3D SPH simulations produced with the MAGMA code (Ross- 
wog & Price 2007), for neutron star components with ini- 
tially no spins, co-rotating spins, and counter-rotating spins. 
The main virtues of this work are 1) the solution of the ra- 
diation transport problem for v e , v e , and "v^" neutrinos for 
eight energy groups, coupled to the hydrodynamical evolu- 
tion of such mergers over a typical timescale of > 100 ms af- 
ter the two components come into contact, and 2) the first 
quantitative assessment of baryon-pollution by the neutrino- 
driven wind produced by such SMNSs. Interestingly, our 
results on neutrino signatures from such merger evens con- 
firm the broad adequacy of previous work that avoided solv- 
ing the radiation transport problem by designing a neutrino- 
trapping, or leakage, scheme (Ruffert et al. 1996, 1997; Ross- 
wog & Liebendorfer 2003). The only noticeable discrep- 
ancy is the much stronger "v^' luminosities we predict with 
VULCAN/2D, something we associate with nucleon-nucleon 
bremsstrahlung not currently accounted for in the above leak- 
age schemes. 

At 10 ms intervals and for each BNS merger model, we se- 
lect a sequence of VULCAN/2D snapshots computed with the 
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FIG. 19. — Variation of the log of the cumulative annihilation rate 
jfmax dyQ+typ) with polar angle 6 raax in the BNS merger models initially 
with no spin (solid), co-rotating spins (dotted), and counter-rotating spins 
(dashed), and shown here at a time of 10 (black) and 60 ms (red) after the 
start of the VULCAN/2D simulation. In other words, S max is the half-opening 
angle of the cone that defines the integration volume for the quantity plotted 
along the ordinate axis. We use the S n results computed with 8 i9-angles and 
accounting only for a density cut for the sites of emission and deposition. 
Energy deposited in the cooling region is not subtracted off, yielding values 
larger by a factor of two compared with the corresponding values given in 
Table 1. 

MGFLD solver and post-process them with the multi-angle S„ 
solver, relaxing the radiation variables, but freezing the fluid 
variables. Based on a knowledge of the energy-dependent, 
species-dependent, and angle-dependent neutrino specific in- 
tensity I v {e v ,n), we compute the neutrino-antineutrino anni- 
hilation rate using a new formalism that incorporates vari- 
ous moments of I v (e v ,n). We find that the total annihilation 
rate computed with the S„ solution and our new formalism 
is larger, but at most by a factor of a few, compared to the 
results based on a combination of a leakage scheme (Ruffert 
et al. 1996) and paired-cell summation (Ruffert et al. 1997). 
With density cuts alone, the annihilation rates based on the S„ 
solution increase by a factor of about two. In our simulations, 
we find that all neutrino luminosities decrease after peak, re- 
sulting in a decrease in the annihilation rate with time, i.e., 
oc f" 1,8 . We find a cumulative total rate that decreases over 
100ms from a few xlO 50 to ~10 49 ergs _1 . Adopting an av- 
erage energy of %kT ~40 MeV for the annihilating v e v e (Bur- 
rows et al. 2006), the number of e~e + pairs produced varies 
over that time span from a few x 10 54 down to ^10 53 s" 1 . 

We show in Fig. 19 the cumulative annihilation rate com- 
puted with the S„ solution and our new formalism at 10 (black) 
and 60 ms (red) after the start of the VULCAN/2D simula- 
tions. Only ~1% of the total is deposited along the rotation 
axis and in a cone with an half-opening angle of ~10°, yield- 
ing a rate of ^10 49 erg s" 1 at 10 ms but down to ~10 erg s" 1 
at 60 ms. Integrating over a tenth of a second, we obtain a 
total energy deposition of <10 48 erg. This result is about an 
order of magnitude smaller than obtained by Setiawan et al. 
(2004, 2006), although in their approach, neutrino emission 
is considerably boosted by shear-heating in the differentially- 
rotating disk. Their annihilation rates are stronger and weaken 
slightly more gradually with time, i.e., oc f _I 5 . In the VUL- 
CAN/20 simulations presented here, we have no physical vis- 
cosity to include shear heating, nor is there adequate resolu- 
tion to simulate the magneto-rotational instability. Our tem- 



peratures and neutrino emissions are therefore underestimates 
of what would occur in Nature. With viscous dissipation as an 
energy source, we anticipate that our results might yield more 
attractive powers for the generation of a relativistic e~e + -pair 
jet, an issue we defer to a future study. 

Baryon loading of relativistic ejecta is a recognized prob- 
lem in the production of a GRB, but this is the first time 
the neutrino-driven wind, the baryon polluter, has been sim- 
ulated dynamically. All three BNS merger simulations we 
performed reveal the birth of a thermally-driven wind through 
neutrino energy deposition in a gain layer at the surface, and 
primarily at high latitudes. The wind blows mostly along 
the polar funnel, contained in a cone with an opening angle 
of ~20° within <500km, but widening at larger distances 
to ^90°. The electron fraction of the associated ejecta is 
close to 0.5, along the pole, but is ~0. 1-0.2 along all lower 
latitudes. Low-T e material at near-equatorial latitudes and 
large distances has a radial velocity below the local escape 
speed, and since local pressure gradients are not negligible, 
it is unclear whether it will eventually escape. Overall, it ap- 
pears that if this pre-black-hole phase lasts for 100 millisec- 
onds, <10~ 4 M Q of "r-process material" will feed the inter- 
stellar medium through this neutrino-driven wind. Moreover, 
such baryonic pollution along the polar direction may affect 
any subsequent relativistic ejecta, although the wind from the 
SMNS may last for only a short time, perhaps 100 ms. Trav- 
elling at a tenth the speed of light, it will reach only out 
to 0.01 light-second, thus much shorter than the duration of 
short GRBs. This initial wind phase cannot provide the sus- 
tained confinement needed for the potential relativistic ejecta. 
Although we focus on the phase prior to black hole forma- 
tion, a neutrino-driven wind may blow after the black hole 
forms too, but this time from the surrounding hot and dense 
torus, and driven by charge-current neutrino absorption at the 
disk surface. Importantly, because of the centrifugal barrier 
felt by particles coming in from sizable orbits in the disk, 
these injected baryons should remain away from the polar 
region, while neutrino-antineutrino annihilation would con- 
tinue to contribute along the rotation axis of the SMNS, but 
now in an essentially baryon-free environment. This pro- 
vides an attractive mechanism for confinement, since neutrino 
emission would yield both the power source for the relativis- 
tic ejecta (caused by neutrino-antineutrino annihilation) and 
the confined neutrino-driven wind (caused by charge-current 
reactions). In the future, we will investigate the neutrino 
signatures that obtain in the context where the compact ob- 
ject is a black hole, and only the surrounding torus mate- 
rial radiates neutrinos. Such a disk wind would offer a nat- 
ural confining ingredient for any relativistic ejecta propelled 
along the rotation axis and powered by annihilation of neutri- 
nos and antineutrinos radiated from the disk, as proposed by 
Mochkovitch et al. (1993). 

The huge amount of free energy of rotation, on the or- 
der of 10 52 erg in the supermassive neutron star and its torus 
disk, together with millisecond orbital periods, suggest that 
the magneto-rotational instability could play a major role in 
redistributing angular momentum, leading to solid-body rota- 
tion, and increasing the magnetic pressure in the correspond- 
ing layers by orders of magnitude. Magneto-rotational effects 
should be strong and would lead to a considerably stronger 
neutrino-driven wind, as in the AIC of white dwarfs (Dessart 
et al. 2007), and delay the formation of a black hole. Simi- 
larly, Dessart et al. (2008) found that in collapsar candidates 
(Woosley & Heger 2006; Yoon et al. 2006; Meynet & Maeder 
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2007), the large amount of free energy of rotation in their iron 
core at the time of collapse can lead to a magnetically-driven 
explosion, associated with very large mass loss rates which 
may compete with the mass accretion rate from the torus 
disk and jeopardize the formation of a black hole. By con- 
trast with such collapsar models, supermassive neutron stars 
formed from merger events are unambiguously endowed with 
a large rotational energy budget (stored in the orbital motion), 
of which a large fraction is differential and may be tapped. 
Hence, the phase prior to black hole formation should be 
modeled at very high resolution and in combination with neu- 
trino transport to address these issues. The maximum neutron 
star mass allowed by the EOS will ultimately determine how 
much mass accretion can take place from the torus-disk, and 
it thus represents a central question for short-duration GRBs. 
Presently, there is a lack of consensus that leaves this issue 
largely unsettled. 
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APPENDIX 

In this Appendix, we broaden the discussion on the VjVj annihilation process presented in § 4, by moving away from BNS 
mergers, and focusing instead on 100-200 ms-old PNSs, formed from the collapse of their Fe or O/Ne/MG core as it exceeded the 
Chandrasekhar mass. Such PNSs are hot because they are young, by contrast with BNS mergers, whose SMNSs are "rejuvenated" 
through shear heating and shocks associated with the coalescence phase. Specifically, we investigate the postbounce neutrino 
heating phase of the 20 M Q progenitor model of Woosley et al. (2002), adopting initally no rotation (model s20.nr) or differential 
rotation (model s20.7r; the initial angular velocity at the center is 7rrads _1 , equivalent to an initial central period of 2 s). These 
models were evolved with the MGFLD and the S„ solvers using the radiation hydrodynamics code VULCAN/2D, and the results 
were described in detail in Ott et al. (2008). To investigate the effect of very fast rotation in a PNS context, we also investigate 
the postbounce neutrino heating phase of the 1.92 M Q AIC model of Dessart et al. (2006b). To make a meaningful comparison 
of the VjVi annihilation rate in all three models, we take their properties at 160 ms after bounce, thus at times when the neutrino 
luminosities are comparable, and, thus, it is the shape of the radiating PNS that differs between models (for completeness, we 
also include some quantities at the last computed times - see Table 2). In this respect, our choice of models extends from the 
quasi-spherical s20.nr model, to the mildly oblate s20.7r model, and finally to the strongly oblate (with a disk-like extension) AIC 
model. In all three models, the S„ neutrino radiation field was computed with 16 i?-angles. Including the ip direction, a total 
number of 144 angles were treated, following the prescription described in Ott et al. (2008). 

In Fig. 20, we present 2D colormaps of the spatial distribution of the energy deposition by v e v e annihilation in the three model 
snapshots considered here. We draw as white the regions of high density (p >10 n g cm" 3 ) where equilibration via charged-current 
interactions is fast, as well as regions in which cooling by charged-current interactions dominates. Any energy input by neutrino 
pair annihilation in these regions is overwhelmed by cooling through charged-current processes and does not contribute to the 
net neutrino gain. Note also that we do not show the distribution for the "z^P^," because it is qualitatively similar to that of v e v e , 
being merely weaker, i.e., with a volume-integrated magnitude that is a factor of about ten smaller. 

In the non-rotating model s20.nr, the Q + {v e v e ) distribution is approximately spherically symmetric, tracing regions of high 
neutrino energy density in the neutrino-diffuse postshock region, and peaking near the lower edge of the gain region around 
^95 km. The v e v e process dominates over that for "v^v^" because the benefit of the favorable "v^' neutrino luminosity is more 
than cancelled by the ^5 times smaller "v^v^" cross sections and by the smaller decoupling radii for "v^" neutrinos, yielding 
systematically smaller collision angles. 

In the 160 ms postbounce snapshot of model s20.nr, Q + {viVi) falls off radially as ~ r" 8 3 (both for v e P e and the "v^v^" anni- 
hilations), which is slightly steeper than the theoretical estimate r~ 8 for r 3> R v made by Goodman et al. (1987). We find and 
list in Table 2 volume-integrated energy deposition rates of 1 .30 x 10 49 and 0.35 x 10 49 erg s" 1 for v e v e and "v^v^" annihilation 
processes, respectively. These rates are more than two orders of magnitude smaller than those associated with the net gain from 
charged-current reactions, which amount to 2.14x 10 51 erg s" 1 at that time. At 500 ms after bounce in the s20.nr model, the energy 
deposition rates have declined to 0.77 x 10 49 and 0. 12 x 10 49 erg s -1 for v e v e and "VyV^ annihilation, respectively. At this time, 
we measure a more rapid radial fall-off oc r" 9 5 for radii <100km and somewhat shallower decline oc r" 8 for larger radii. These 
results echo the steep decrease with time of the annihilation rate found for the BNS mergers discussed in this paper, suggesting 
that neutrino-cooled compact objects are quickly weakening neutrino-annihilation power sources in the absence of energy sources 
such as shear heating in a differentially-rotating disk (Setiawan et al. 2004, 2006). 

Our results for model s20.nr support the conclusions of previous studies (Cooperstein et al. 1986,1987; Janka 1991) that argued 
that ViVj annihilation contributes little to the neutrino heating in quasi-spherical postbounce supernova cores. GR effects (bending 
of neutrino geodesies and redshift), which yield larger annihilation rates for very compact configurations (Salmonson & Wilson 
1999; Bhattacharyya et al. 2007; Birkl et al. 2007) are most likely irrelevant at the large radii at which annihilation may have any 
significance in our models. GM/Rc 2 is already as low as ~0.02 at the inner edge of the polar gain region in models s20.nr and 
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FIG. 20. — Colormaps of the energy deposition rate Q + (u e v c ) by electron neutrino pair annihilation in the nonrotating model s20.nr (left), the rotating model 
s20.7r (center), and in the 1.92-Mq accretion-induced collapse (AIC) model of Dessart et al. (2006b) (right), and in each case at 160ms after bounce. The 
energy deposition by the "v^v^" annihilation process exhibits the same qualitative spatial distribution and is not shown here (the corresponding cumulative 
rates are typically a factor of ten weaker than for v e v e ). The black lines demarcate regions in which charged-current losses dominate and regions deep inside 
the PNS where /3-equilibrium prevails. Any energy deposition from pair annihilations are instantly re-emitted in these regions by charged-current interactions 
and no net energy deposition results. In the snapshots shown here, we find integral net energy deposition rates by neutrino pair annihilation of 1.65xl0 49 , 
1.44x 10 4 ' and 0.59x 10 45 erg s -1 for models s20.nr, s20.7r, and the AIC model, respectively. In this computation, we exclude the regions with a density larger 
than 10" gem -3 , as well as regions of negative net gain (shown in white). These numbers are one to two orders of magnitude smaller than the corresponding 
rates due to charged-current neutrino absorption. Note that all models are run with 16 i9-angles. 



TABLE 2 

Cumulative ^,p,- Annihilation Rates 



Model Name 


t-tz, 


£(cc) 


E(v e V e ) 




£('VpiV) 




Polar 




(ms) 


(10 49 ergs s-') 


(10 4s ergs s 


-') 


(10 49 ergs s" 1 ) 


Fall-Off 


s20.nr 


160 


214.0 




1.30 


0.35 


oo 




s20.nr 


500 


81.8 




0.77 


0.12 


oc 


w .-9.5_ r -8.0 


S20.7T 


160 


160.2 




1.03 


0.41 


00 




S20.7T 


550 


26.2 




0.29 


0.12 


oc 




AIC 1.92-M 


160 


53.7 




0.42 


0.17 


OC' 




AIC 1.92-MfT) 


773 


48.9 




0.19 


0.02 


oc 





NOTE. — Summary of the cumulative (but instantaneous) neutrino-antineutrino annihilation rate calcula- 
tions, t-tfc corresponds to the post-bounce time of the computation. £(cc) is the integrated gain from charged- 
current interactions, while E(y e v e } and E^'v^P^") denote the integral energy deposition by annihilation of 
v e v e and "u^u^" pairs (accounting for v^v^ and v T u T ), respectively. The rightmost column gives the ap- 
proximate power-law exponent that describes the radial fall-off of the energy deposition rates along the polar 
direction (we see the same distribution for both the v e v e and "v^v^ cases). For a spherically-symmetric 
model in which neutrinos irradiate isotropically from a neutrinosphere, Q + (Vii'j) oc r~~ s at radii large compared 
to the neutrinosphere radius (Goodman et al. 1987). The ebbing annihilation rate with time results from its 
strong J J' dependence (eq. 11), and reflects PNS cooling. [See text for discussion.] 

s20.7r and not larger than ~0.06 in the AIC snapshot. 

Model s20.7r is rapidly rotating and has large pole-equator neutrino flux asymmetries (Ott et al. 2008). Though prolate and 
quickly varying with radius, the <2 + (iw) distributions in model s20.7r do not vary with angle by more than a factor of 2 (neglecting 
for now the large equatorial cooling regions between ^100 and 200 km). Along the poles, charged-current losses negate net 
energy deposition by viv-, annihilation at radii <80km. Again, Q, + (v e v e ) dominates over Q^i^v^v^"), while both decrease with 
radius by ~ r" 7 8 at 160 ms after bounce and ~ r" 6 5 at 550 ms after bounce. As listed in Table 2, we find at 160 ms after bounce 
volume-integrated energy depositions of 1.03 xlO 49 and 0.41 x 10 49 erg s -1 , for v e v e and 'V^i^," respectively. At 550ms after 
bounce, these values have decreased to 0.29x 10 49 {v e v e ) and 0.12x 10 49 ergs s" 1 ('V„i>„"). The total annihilation contribution to 
neutrino heating is equivalent to ^1% (~1.5%) at 160 ms (550 ms). 

The right panel in Fig. 20 shows the energy deposition by v e v e annihilation in the AIC model at 160 ms after bounce (see Dessart 
et al. 2006b for details). In this rapidly rotating oblate PNS (the central period is ~2.2 ms at this post-bounce time), neutrinos 
are emitted from both the central object and its extended moderately hot equatorial disk-like structure (more specifically from 
the pole-facing side of these side lobes). Significant energy deposition by neutrino pairs occurs at small radii in a wide polar 
wedge of ^60° (widening with radius to ~120°) and Q + (v e v e ) reaches peak rates near the PNS surface at ^15 km of up to 
10 29 erg cm 3 s with <2 + (" l V I V") being globally smaller by an order of magnitude. The radial decline of £? + (z-w) along the 
polar direction becomes shallower with time, with power-law exponents varying from ~-8.3 to ~-5.8 from 160 to 773 ms after 
bounce. 

In this AIC model, both Q + {v e v e ) and Q^i^v^v^') reach local values at small radii that are of the same order of magnitude as 
the peak values in energy deposition per unit volume by charged-current interactions, yet the very limited volume of the high- 
Q + (i>iVj) gain regions leads to only very modest integral values at 160ms after bounce of 0.42 xlO 49 and 0.17xl0 49 ergs _1 for 
v e v e and "v^v^' pair annihilation, respectively, in our Newtonian model. At 773 ms after bounce and in the same order, these 
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values are 0.19x10 and 0.02x10 ergs . Compared with the charged-current interactions driving the post-explosion wind 
phase of the AIC (Dessart et al. 2006b), even a ten times larger integral E{v;i>i) (e.g., via GR effects) would still amount to only 
< 10% of the total charged-current energy deposition rate in this model. 
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